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FOREWORD 


The program described herein was conducted by the General Electric 
Company Aircraft Engine Business Group under NASA Contract NAS3-23525. 

Mr. Robert R. Tacina of the Aerothermodynamics and Fuels Division, NASA Lewis 
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Investigator; D.W. Bahr, Technical Program Manager; E.J. Rogala, Program Man- 
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ment of the axisymmetric program GETREF; D.L. Burrus directed the experimental 
program, and with the assistance of P.H. Waymeyer, performed the remaining 
assessments; and W.C. Colley assisted in defining needed computer programming 
changes in the 3-D elliptic model. 

The Final Report is divided into two volumes, with Volume II containing 
the experimental data obtained under this contract. An Executive Summary 
report was also prepared. 
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1.0 SUMMARY 




The Aerotherinal Modeling Program - Phase I was conducted as part of the 
overall Hot Section Technology (HOST) Program, initiated by NASA. The Phase I 
program was a 9 month effort involving the assembly of a benchmark quality data 
b^se from selected available literature and from General Electric's engine and 
combustor component test data. Additional definitive data were obtained from 
an experimental test program conducted as part of the Phase I effort. This 
experimental effort involved the measurement of temperature distributions for 
a series of combustor- like flows of progressively increasing flow complexity. 

Using selections from this assembled data base, assessment studies were 
cojnducted to determine the predictive accuracy of, and identify deficiencies 
within, the various modules comprising the overall aerothermal model used at 
General Electric. A major need identified was for an internal flow module for 
improved hot gas side inputs to the heat transfer module. Assessment of the 
internal flow module was conducted with 2-D parabolic and elliptic, as well 
as I 3-D elliptic internal flow calculations of definitive test data selected 
from the assembled data base. Calculations of combustor liner metal tempera- 
tups, pessure loss performance, and airflow distribution were performed 
using the aerothermal modules in general use at General Electric. These cal- 
culations wpe compapd with test data selected from the data base to assess 
the predictive capabilities of these modules and identify their deficiencies. 

The 2-D axisymmetric assessments provided a methodical examination of 
the mathematical techniques and the physical submodels, while the 3-D assess- 
ments focused on overall usability as a design tool. 

The axis 3 rmmetric assesments of the mathematical/physical submodels indi- 
cated the following: 

• For present purposes; the probability density function, pdf, is a 
satisfactory treatment of the turbulence/ chemistry interaction. 

• Accurate predictive capabilities require improved accuracy in the 
numerics by way of an alternate to central-differencing for 3-D com- 
bustor problems. 

• The K-e prbulence model is imperfect, but computationally feasible 

alternatives, except for swirl or streamline curvature corrections 
are not available. * 

• Calculated results are sensitive to the boundary conditions which 
need to be better defined. 

The assessment of the usability of the available 3-D codes showed: 

• A variety of difficulties in utilizing the available codes. 

• High computer time running costs. 






Some juggling of input parameters, not necessarily consistent with 
the physics, was often required to match the data. 

The two available codes showed similar penetration and subsequent 

ferenf ^ input parameters, although dif- 

field^ calculated for the turbulence throughout the flow 

temperature pattern from a com- 

Ho!evL ? approximately reproduced by the calculations. 

?mnro« M methods of defining the initial dome flows and 

improved numerics are needed before the calculations are sufficient^ 
accurate to be a useful design tool for durability analyses ^ 
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2.0 INTRODUCTION 
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ature le^el with^?L«i ! ® ^gn^ficant increase in turbine' inlet temper- 
u ^ attendant thermodynamic cycle benefits. Further cycle 
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hill,.*,, f -w acnieved with higher pressure ratio compressors The 

higher turbine inlet temperatures and compressor pressure ratiLwit^^nr: 
re spending higher temperature cooling air has created a hostile environment 

maL^eL^° components. As a result, a major share of the engine 

maintenance cost is associated with the combustor and the turbine airfoils. 

costs^°NTsrh« «®=tion maintenance 

costs, NASA has initiated the Hot Section Technoloev n i 

element of this overall program is the Hot S^Jt^Jn^L EogrAlrotherSa! " 

vaJidaS objective of this program i^L eJowHud 

validate improved analysis methods for use in the design of aircraft turbinA 

exnerLdT ''°“; combustor analysis capabilltiL can be 

teristico significant improvement in the life and durability charac- 

a sienificant red turbine components. Another benefit Luld be 

mentf? “^“‘^tion in expensive full-scale development testing require- 


etaiiea internal flow module under de/elopment. 
of th^PhLe program was initiated in 1982. The objectives 
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accuracy of these calculations was determined through comparisons with the 
data base. Based on this assessment, specific deficiencies in the generic 
model were identified, and recommendations for improvement were formulated 


3.0 PROGRAM SCOPE 


3.1 OBJECTIVE 

The objective of Phase I of the Aerothermal Modeling Program was to iden- 
tify deficiencies in the aerothermal model that is currently used or planned 
for use at General Electric in the design of gas turbine engine combustors. 

This effort was accomplished by assessing the predictions of exising modules 
with benchmark quality test data. This effort, together with other similar 
efforts sponsored by NASA, was intended to assess the current state of the art. 
It was also intended to define Phase II efforts needed to evolve improved 
and more advanced aerothermal analysis methods, thereby enhancing available 
capabilities for the design of high performance and durable combustors for 
advanced aircraft turbine engines. 

j This Phase I effort was focused on the ability of the existing model to 
predict or analyze those features associated with durability. These primarily 
are: (1) the combustor liner metal temperature distribution, and (2) the 

cbmbustor exit gas temperature patterns, which in turn affects the life of the 
downstream turbine components. Other features of interest to the combustor 
designers, such as blowout limits and pollutant emission levels, do not 
directly affect design for durability. 

3.2 THE GENERAL ELECTRIC AEROTHERMAL MODEL 

Figure 1 shows schematically General Electric's calculation modules within 
the overall Aerothermal Model. Except for the internal flow module, the mod- 
ules are well developed and have been in use for many years. The internal flow 
module is planned for use in the design process at General Electric after it 
has been developed and/or demonstrated to have useful accuracy. Thus, assess- 
ment of the candidate internal flow modules available at General Electric con- 
stitute a major portion of the efforts reported herein. 

The model is described in detail in Section 5.0 of this report. Briefly 
the modules shown in Figure 1 provide the following calculated results. The 
Preliminary Design Module provides the basic features of the combustor, 
including its flowpath and the number of fuel nozzles, early in the design 
process. The Diffuser Module calculates diffuser pressure losses and the 
three individual pressure levels that feed the inner and outer flowpaths and 
the dome region. The Flow Distribution Module, using the diffuser exit 
pressure levels, total airflow, fuel flow, and the combustor aperture areas, 
ahd features, calculates the flow and flow angle through each aperture, as 
well as a one-dimensional internal gas temperature and the overall pres- 
sure drop. The program has features to determine effective flow areas from 
the particular aperture shape and approach flow velocity. A large data base 
exists at General Electric for providing effective area correlations and new 
measurements are made on individual apertures whenever needed. 


The Heat Transfer Module utilizes the flow distribution, including liner 
cooling slot flows, the velocities on the cold side of the liner, and some 
modification of the one-dimensional internal gas temperatures and velocities 
from the Flow Distribution Module. Through conventional heat transfer equa- 
tions, together with correlations of cooling film effectiveness based on wind 
tunnel data, the metal temperatures are calculated throughout the combustor 
metal structure. Figure 2 illustrates the type of node pattern for which tem- 
peratures might be provided for one of the several cooling panels along the 
length of a combustor liner. This information, together with aerodynamic 
pressure loading from the Flow Distribution Module, is then used as input to 
stress and life analysis procedures to provide life estimates or to indicate 
the life effects of a postulated design change. The input also includes 
temperature distribution information in the circumferential direction. 

Metal temperatures are measured on combustors, but in only a finite number 
of places. The available stress and life analysis procedures need more than 
isolated temperature measurements. They need the detailed temperature distri- 
bution provided by the Heat Transfer Module. The major source of error in the 
Heat Transfer Module is in the estimation of the inputs on the hot gas side; 
the local flow velocities, gas temperatures, and flame radiation levels. 

If an accurate 3-D elliptic computer code for calculating the internal flow 
field existed, it could provide substantially better inputs. 

At General Electric, the combustor exit gas temperature patterns are 
measured in detail, typically at seven radial positions and every 1-1/2 
degrees around the circumference. This information is of adequate detail for 
the turbine designers . However, improving or modifying this temperature 
pattern during the development of the combustor usually involves many cut-and- 
try test efforts. Changes are often defined with the help of existing data or 
correlations of simple jet penetration experiments in a cross flow oi measure- 
ments of flows emerging from the swirl cups surrounding a fuel nozzle in a 
simple bench experiment. An accurate internal flow module would however be 
very valuable in providing more detailed and precise inputs and, thus, in 
reducing total development tests. 

Four basic codes for calculating elliptic reacting flows are available at 
General Electric: (1) the code assembled by the Garrett Turbine Company under 

the Army-sponsored Combustor Design Validation Program efforts, Mongia et al 
I979I*. ( 2 ) the very similar code assembled by the Northern Research and 
Engineering Corporation, NREC Report No. 1420^; (3) the TEACH code developed 
initially at Imperial College which uses essentially the same methods as the 
previous two codes but in a two-dimensional instead of three-dimensional frame- 
work; and 4) the two-dimensional (planar or axisymmetric) research code, GETREF 
(General E^lectric Turbulent Reactant Flow) developed at General Electric, A 
parabolic version of this model has also been developed for use in shear flows. 

Prior to this effort, plotting routines for the output and some improved 
flexibility for input had been coupled at General Electric with the three- 
dimensional elliptic codes and this new format containing both codes was 
defined as the INTFLOW program. Creation of new programming was not to be 
part of the scope of effort on this Phase I Aerothermal Modeling program. 

Superscript numbers are reference numbers. 
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The two-dimensional GETREF code provided a valuable addition to the three- 
dimensional codes because the physical submodels and numerical techniques are 
more rigorous. Figure 3 is an illustration of a calculation from the GETREF 
code dealing with the initial swirl cup portion of a combustor where the flows 
are initially axisymmetric. For a typical full annular combustor, a model of 
a segment containing one fuel nozzle is modeled. Figure 4 illustrates a 
calculated flow field. The length of the arrows are proportional to the 
calculated velocities. Only recently have computer storage capabilities 
been able to accept a model with this much grid detail, see Figure 5. Little 
additional grid detail can be introduced at this time due to computer storage 
limits, but the cost of obtaining a converged solution with more grid could 
become prohibitive at present computer running time costs. 

Some of the immediately obvious limitations of the existing 3-D internal 
flow codes that were recognized at the outset and influenced the assessments 
conducted in Phase I were, as follows: 

1. With the grid shown in Figure 5 the so called "hybrid" numerics in 
the program involve the introduction of significant artificial or 
numerical diffusion. This produces mixing rates greater than speci- 
fied by the turbulence model. Whenever the individual calculation 
cells have a Peclet number above 2.0, second order accurate central 
differencing numerics cannot be used and upwind differencing (first 
order accurate) is used which introduces significant numerical diffu- 
sion. Figure 6 shows that the cell Peclet numbers are below 2.0 in 
only isolated regions. Thus, without the ability to use much finer 
grids, the model inherently involves significant numerical diffusion 
for combustor flows. Claus, 19833, showed that even with much more 
detailed grid for a single dilution hole than possible for a real 
combustor, analysis inaccuracy results from the numerical diffusion. 

2. In order to introduce noncyclindrical combustor walls, the stair- 
step technique must be used and this involves grid concentration in 
unneeded regions and crudeness in calculating film-cooled walls. 

3. In the existing codes, the geometry-based swirl cup and fuel injec- 
tion inputs do not correctly produce the initial fuel distribution 
for all combustors. Swirl cups have been tested at General Electric 
with minor geometry changes- that make major differences in metal hot 
streaks and exit gas temperature patterns; these geometry changes 
would not result in any change in the model input unless some 
knowledge (external to the computer program) of the nature of the 
flows emerging from the swirl cups was available. 

Accordingly, improvements addressing these deficiencies must be evolved 
before the codes become the effective design and analysis tools suitable for 
defining combustors with enhanced durability. 
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Figure 3, Velocity Field Calculation Using GETREF Code (2-D). 
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Figure A. Calculated Velocity Field for a Combustor Using 3-D Elliptic Code (K Plane 15). 
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Figure 6. Calculated Peclet Number In the X-Dlrectlon from 3-D Combustor Calculation 






3 . 3 DATA SELECTED FOR COMPARISON WITH MODELS 


Data for assessing the accuracy or deficiencies of the modules in this 
effort are available from three sources: (1) a General Electric aircraft tur- 

bine engine combustor, (2) experiments reported in the literature, and (3) 
experiments conducted at General Electric in this effort. The engine combus- 
tor data include measured metal temperatures and exit gas temperatures, as 
well as measurements of the flows and pressures through the diffuser and the 
annular passages surrounding the combustor. While the temperature data may 
be sufficient for a check of the overall validity of the modules, these data 
alone do not provide sufficient detail to permit the individual sources of 
error within a given model to be identified, particularly the internal flow 
codes. 

On the other hand, geometrically simpler experiments that contain the 
same physical phenomena permit detailed probing of the flow field. Data 
obtained in such experiments are, therefore, defined as a "benchmark" quality 
data for use in checking the mathematical physical elements of the modules. 
Most of these experiments are two-dimensional for ease of data taking. This 
permits a more accurate assessment of the physics since, for example, numeri- 
cal inaccuracies can be minimized by using more grid points per coordinate 
direction. 


Data from four axisymmetric "benchmark" experiments were chosen to assess 
the mathematical physical model. The first is an existing experiment performed 
at General Electric and consists of an axisy^etric fuel jet burning in a 
coflowing air (Lapp, et al., 1982^). The primary reasons for its selection 
are: (1) the removal of "numerical diffusion" errors since the governing 

equations become parabolic and; (2) the availability of high-quality laser 
doppler velocimeter LDV and spontaneous Raman scattering data. An 
unambiguous evaluation of the turbulence-chemistry interaction model can, 
therefore, be made. In addition to the jet experiment, three axisymmetric 
recirculating flow experiments were chosen: an isothermal annular dilution 

jet experiment (Green and Whitelaw 19805), an isothermal double-swirled pipe 
flow (Vu and Golden 1980^), and a bluff-body stabilized non-premixed flame 
(Liightman, et al.^). These provide checks of the model under progressively 
more complex circumstances. The above four experiments were analytically sim- 
ulated using codes (elliptic GETREF or its parabolic version, as relevant) 
developed at General Electric. These codes are described in more detail in 
Sections 6.0 and 7.1. 


For simple round dilution jet flows in a crossflow, a three-dimensional 
flbwfield, data from Walker and Kors 1973,® was chosen for examination. In 
the past, these data provided a large matrix of penetration data, permitting 
useful penetration correlations. Hoi deman and Walker 1977°. These data and 
correlations have been used in the past at General Electric in combustor 
design and development work. 


To provide data from more complex flowfields, a series of experiments was 
devised that involved progressively more complex flows beginning with a row of 
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jets in a crossflow and followed by opposing jets, alternate jets of different 
sizes, and a second row of dilution jets with the two row pattern simulating 
the dilution pattern of one General Electric combustor. These initial tests 
in the series were conducted with cylindrical inner and outer walls and a flat 
uniform inlet condition. The two available three-dimensional elliptic computer 
codes use cylindrical coordinates and, as received at General Electric, did 
not have capabilities for inputting combustor walls with a converging exit 
area as is present in actual General Electric combustors. This input capabil- 
ity to the INTFLOW program has been introduced at General Electric along with 
the ability to input swirler flow in an annular combustor. However, the 
experiments were devised to progress to provide data, without the complexity 
^f noncyl indr ical walls and with simple uniform flow simulating the combustor 
dome flow, first with one row of dilution holes, then with two rows simulating 
the FlOl dilution pattern. Tests with recirculation at the dome location and 
the two rows of dilution jets were introduced in two steps. First, a nonswirl- 
ing dome flow was used with flow axisymmetrically diverging from the fuel noz- 
zle position followed by a swirling flow. This last case incorporated a dome 
with counterrotating swirl cups like those used in a specific General Electric 
combustor. Then the configuration was tested with the true wall and dome 
contours of this combustor. All of the above tests were conducted without 
injection. As in the Walker and Kors® data for the simple dilution jet, 
the approach How was introduced at a different temperature than the dilution 
jets and detailed temperature measurements were the primary data. The experi- 
ments were then extended to include combustion, first with a gaseous fuel that 
avoids consideration of fuel drop atomization and vaporization, and then with 
a liquid fuel. The effect of film slots was introduced in these last tests. 

The data was supplemented with nonburning data taken on the swirl cup and fuel 
nozzle by itself. 


3.4 PROGRAM APPROACH 

The approach for assessing the modules already in general use at General 
Electric was accompished by relatively straightforward comparisons with avail- 
able General Electric data. 

The approach for assessing the Internal Flow Module was more complex and 
is discussed in the paragraphs below. 

The assessment of the usability of the existing internal flow codes was 
conducted in two parts. The first part consisted of comparing "benchmark" 
quality data with the predictions from the axisymmetric elliptic research code 
(GETREF) and its parabolic version as applicable. These comparisons permit 
very methodical examinations of the mathematical techniques and the physical 
submodules. Based on the results, general conclusions regarding future module 
development were then drawn. 

I The second aspect of the assessment addressed the usability of the 3-D 
codes. In these assessments, the methodical assessment of numerical error and 
the physical submodules as conducted with the axisymmetric code GETREF was not 
attempted. The assessment approach was more the approach of the eventual 


engineer user rather than the original progranuner. The objectives, therefore, 
included observation of the ease of input, convergence rate, computing time, 
and of the accuracy of the results for a typical user. Primary focus was 
placed on the ability of the codes in two areas: (1) exit gas temperature 

prediction and (2) hot gas side inputs to the Heat Transfer Module. 

For exit gas temperatures, it was postulated at the outset that, even if 
the primary combustion zone flow fields were poorly predicted a useful tool 
could still be develoepd if the dilution jet mixing process could be ade- 
quately calculated. An approximately correct dome flow analyses could be 
input and the subsequent dilution mixing calculation could be used to find the 
relative effect of some postulated dilution hole change. 

‘ In the case of the hot side inputs to the heat transfer calculation, 
again design usefulness could be expected even with somewhat limited accuracy 
of the internal flow calculation. First, the protection of the wall by the 
film cooling could still be dealt with by the present methods which involve 
the correlation of wind tunnel data for many film slots. Hence, only the hot 
gas side data just beyond the boundary of the film protection layer would be 
needed. It was recognized at the outset that the existing codes did not have 
the capability to predict flame radiation levels better than current estimating 
methods. The development of this capability and, hence, its evaluation was 
relegated to some future effort. Again, as in the case of the gas exit tem- 
perature patterns, an improved technique for inputting the air swirler and 
fuel flows into the dome would be needed to obtain this capability. 

For real combustors, the basic flow field may be controlled much more 
by the pressure fields than by the turbulent shear layers prominent in simpler 
axisymmetric experiments in the literature. For example, the length of a 
recirculation zone is a critical parameter in the axisymmetric experiments. 

In a combustor, the recirculation zone length frequently ends just upstream of 
the first dilution hole location. Its length is predominantly controlled by 
the pressure field created by dilution holes rather than the turbulent mixing 
phenomena in the axisymmetric experiments. Thus, it may be possible to 
calculate reasonable combustor flow fields even though turbulent mixing is not 
being computed accurately; turbulence models with improved accuracy may not be 
essential. 

In the effort herein it was not presumed that combustor models needed to 
be as accurate or as detailed in all aspects as was needed for good compari- 
sons with the simplified experiments. Thus, the efforts with the 3-D models 
weire to see if useful results for durability design and development could be 
achieved even though it was recognized that the node detail available would 
be inadequate by the standards explored for the axisymmetric experiments. The 
number of nodes utilized was up to the limits available on General Electric's 
computer, an IBM 3081D. 


4.0 STRUCTURE OF THE REPORT 


Sec tion 5.1 describes General Electric's aerothermal model. 

data availabla in tha 

JLiceracure from which four experiments were selected to helo assess thp mrtHpi 
Comments on modeling from the literature are also included. In addition 
discussion of engine and combustor component test data available at General 
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describes the test program conducted at General Electric 
under this effort to provide additional data. ‘general Electric 

Gener arklecrri-c the approach in modeling internal flow field at 

General Electric. The discussions include simple test examples on numeric 

schemes that are possible alternates to the current hybrid Lheme of central 
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S^ tion 7 2 compares results of the 3-D model assessment with benchmark 
experiments selected from the literature, and with experimental test data 
obtained as part of the Phase I program. 

Se ction 7.3 compares results of the assessment of the aerothermal mod.,1p« 
in general use at General Electric with selected General Electric test data! 

Sec tion 7.4 summarizes the assessments of the aerothermal model. 

discusses the overall conclusions obtained from the various 
assessment studies conducted. Recommendations for improving current aero- 
thermal models are developed and presented in the form of a prioritized list. 

Referencejui^ contains the list of , references used in this report. 
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5.1.1 Preliminary Design Module 
Description 

! The combustion system preliminary design module (COMBIO), calculates 
^he overall dimensions and preliminary contours of the combustion system 
flowpath design. This computer program is used to generate the dimensions 
and coordinates required to calculate the key combustor design parameters 
and combustor performance estimates (Figures 8 and 9). This preliminary 
design estimate serves as a starting point for the more detailed design 
and analysis modules that are used to finalize the combustion system design. 


Analytical Methodology 

The calculations within the program are based on empirical data correla- 
tions that have been in existence at General Electric for many years. In the 
program values for reference velocity and space rate are varied to provide 
calculations for an array of design cases. A flowchart of COMBIO, the pre- 
liminary design module, is shown in Figure 10. 


5.1.2 Inlet Diffuser Module 
Description 

I The inlet diffuser module is comprised of several unlinked computer pro- 
grams that perform three basic calculations. Prior to entry into the module 
the preliminary overall sizing of the diffuser has been determined in the 
preliminary design module. This information is used along with the aerothermal 
requirements of the overall combustor design as input into the inlet diffuser 
module. 

The initial calculation step involves generating the diffuser wall 
contours, and is performed by the CONTOUR program. The next step in the cal- 
culation procedure involves determining the streamlime distribution through 
the diffuser flowpath and potential flow calculations along the streamlines. 
Computer programs CAFD or STC are used to perform these calculations. The 
velocity distributions along the diffuser wall contours are obtained as a 
result of the potential flow calculations. The velocity distributions are 
used to represent equivalent straight diffusers that in turn are compared 
wijth the GE diffuser flow regime criteria to determine the flow separation 
stjall margin of the design. In the final calculation step in the inlet 
diffuser module, the DIFLO program calculates predictions of the diffuser 
aerodynamic performance. A flowchart for the Inlet Diffuser Module is 
provided in Figure 11. 
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Figure 10. Plow Diagram of Preliminary Design Module 
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Figure 11. Flow Diagram of Inlet Diffuser Design Module. 







Analytical Methodology 


The CONTOUR program calculates the coordinates of the centerline of 
the passage or passages (in the case of multiduct diffuser designs), and, 
for a selected passage area distributions, it calculates the X and Y coordi- 
nates for the passage walls. For a specified number of calculation inter- 
yals, this program generates a table of passage centerline coordinate^, the 
distance along the centerline, the centerline area distribution, and outer 
and inner wall coordinates. The program can handle raultipassage diffuser 
designs. 

The passage area distribution is based on a cubic equation with a one 
parameter shape factor which can vary from 0.500 for an almost linear distri- 
bution, to 0.667 for a "bell-shaped" distribution. 

These contours can be used as input for the potential flow analysis 
calculation step with precise definition of the boundaries, which greatly 
improves the accuracy of the analysis. Also, the output from this program 
can be used to specify the diffuser wall coordinates for the mechanical 
design of the diffuser casings. 

The CAFD program utilizes the streamline curvature solution technique 
tp compute the circumferentially averaged flow properties. Thermodynamic 
relations are based on variable gas properties and are consistent with other 
General Electric data reduction and cycle analysis programs. Axisymmetric 
force fields and thickness blockage effects are used to represent the 
effects of struts and vanes. The annulus shape, circumferential vane or 
strut thickness distribution, and angle distribution are specified. Calcu- 
lation stations can be arbitrarily shaped. Many options are available for 
the user to specify the vane row turning and loss. These include design 
options such as angular momentum and loss coefficient, and analysis type 
options such as flow angle and loss coefficient. 

The cross stream momentum equation is used to calculate velocity gradient 
as a function of (1) total enthalpy and entropy distribution as determined 
from "user" input and application of the energy equation, (2) blade force 
effects determined from blade geometry and work distribution; and (3) local 
slopes and curvature of the streamlines. The slope and curvature are deter- 
mined by a curve fit of the assumed streamline positions. This curve fit 
may be specified as a rotated-spline, or a combination of a spline and a 
three-point parabola. 

An initial set of streamline positions are curve fit to find slopes and 
curvatures, the cross stream momentum equation and the continuity equation 
are applied at each calculation station to determine the shift in the stream- 
line positions. These indicated shifts are based on fixed values of slopes 
and curvatures. Moving the streamlines by the indicated shifts would result 
in gross overcorrection and the iteration would diverge. Therefore, the 
actual streamline movements are determined by a second order correction 
equation. The driving force for this equation is the indicated shifts found 
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in the flow balance. The actual streamline adjustments, as calculated by 
the correction equation, are applied without damping. The new positions are 
used to repeat the above loop. Convergence occurs when the maximum indicated 
difference is acceptably small. The convergence rate of this system is found 
to be quite rapid. 

The major limitation of the CAFD program, as applied to combustor inlet 
diffusers, is that constant pressure free streamlines cannot be calculated. 
Thus, free streamlines downstream of the dumping point of step type diffusers 
are not accurately simulated. 

The STC program (Ferguson and Keith 1975^®) is a very comprehensive com” 
pressible, inviscid, potential flow analysis program developed under NASA 
sponsorship. Within the program an iterative technique is used to solve 
planar or axisymmetric flow fields. The method utilized to compute the flow 
field is the Streamtube Curvature Relaxation technique. 

The STC approach utilizes a number of confluent streamtubes, with slightly 
different properties that are "dded together to obtain the total flow in the 
channel. Each streamtube is handled in much the same way as is the one 
streamtube in the one-dimensional problem. In the limit, as the size of the 
individual streamtubes approaches zero, the STC method satisfies the inviscid 
equations of motion exactly. 

The solution method is an extension of the streamline curvature method, 
and is briefly described as follows: First, a crude grid of streamlines 

and orthogonal lines is assumed. Second, the curvature of the streamlines 
at each of the grid points is evaluated. Third, the momentum equation is 
integrated along a line normal to the streamlines to obtain velocity, and 
the continuity equation is integrated to determine the "correct" streamline 
positions (for the assumed curvature field). Fourth, an adjustment (An) is 
computed by considering: (1) the difference between the computed and assumed 

streamline positions and (2) the effect of the implied curvature modification 
in the integrated momentum equation. Finally, the streamlines are reposi- 
tioned by the An values. 

Because the movement of any one grid point alters the velocity at nearby 
points through a change in curvature, it is highly desirable to account for 
these interrelating point adjustments simultaneously. The utilization of a 
simultaneous solution procedure, employed here, is not part of the classical 
streamline curvature method. 

The streamline curvature method is extremely powerful as indicated by 
the following features: 

• No additional complexities arise when the flow is rotational. 

• The slip line between an exhaust jet and the external flow can be 
handled precisely. (The procedure is to consider two coincident 
streamlines. Their position and pressure are the same; their 
velocity and stagnation properties may be different.) 


• The streamline/orthogonal line grid provides a mapping of the flow 
field into a rectangular domain. ThiS' is helpful from the stand- 
point of computer program organization. 

The STC program has been designed to; 

• Handle multiple streams, 

• Place grid points at locations in the flow field where they are 
needed, as determined by local variations of the dependent vari- 
ables . 

• Allow external flow analysis by incorporating matched near field 
and far field solutions. The far field solutions are obtained 
analytically utilizing small perturbation theory. 

The operations performed by the STC program many be outlined as follows: 

1. Define the flow regions and locate (approximately) the "primary" 
orthogonals and the streamlines which divide the internal and 
external flows. 

2. Refine the grid as required by inserting additional streamlines 
and orthogonal lines between those already existing. 

3. Compute the streamline angle and curvatures. 

4. Compute the orthogonal line angles and move the grid points along 
the streamlines to obtain orthogonality. 

5. Compute the velocities on the "far field" boundary. 

6. Adjust the flow rates in the exhaust streams, if any, to meet the 
calculated choking flow rate. 

7. Integrate along each orthogonal the momentum and continuity 
equations . 

8. Determine if the streamline positions are within a "rough tolerance". 
If so, return to Step 2 for additional grid refinement (unless grid 
refinement limits have already been reached). Otherwise, continue 

to Step 9. 

9. Determine if the streamline positions are within final tolerance. 

If so, jump to Step 13. Otherwise continue to Step 10. 

10. Set up the matrix equation for the streamline correction, 4n. 

11. Solve the matrix equation. 

12. Modify the streamline positions by 4n, and return to Step 3. 


13. Calculate and print the output quantities; then return to Step 1 
for the next case, if any. 

The first operation includes reading the input for a description of the 
gepmetry and flow properties. The first step in the programmed logic is to 
develop a table for orthogonals or calculation stations for each of the 
several flow regions. The regions are determined so the calculation can 
proceed from upstream to downstream. The boundary of each region is defined 
as a primary orthogonal. The initial grid developed contains only the 
primary orthogonals and the double streamlines which separate the various 
streams. 

i The iterative sequence is to start with a crude grid, as noted above, 
and to repeat Steps 3 through 12 until the flow balance error is small. This 
is often accomplished in one or two iterations. The grid is then refined to 
the next level, and the field reconverged. The refinement/convergence pro- 
cess is continued until the grid refinement criteria is satisfied, or alter- 
nately, until storage limits are reached. Greater detail into the analytical 
methodology used in the STC program can be found in Reference 10. 


The DIFLO program calculates total pressure losses and static pressure 
rise coefficients in each passage of the combustor inlet diffuser. An over- 
all, mass-weighted pressure loss for the diffuser system is also determined. 
These losses include prediffuser momentum losses, based on General Electric 
straight annular diffuser test data; flow dumping losses, based on General 
Electric bench-side test results for simular step configurations; and para- 
sitic drag losses of struts, mounting supports, and fuel nozzle systems. The 
losses in the prediffuser section are adjusted to account for velocity profile 
effects. Momentum mixing losses along free streamline boundaries are also 
included. 

' Predictions of combustor inlet diffuser pressure losses are based on 
dilffuser test data. A large number of full-scale, full annular diffuser 
test programs have been conducted by General Electric over the past 25 years 
for a variety of different combustor inlet diffuser configurations. The 
extensive data base accumulated^ from these test programs and from other 
sources has been used to develop the DIFLO computer program. 


Predictive Capabilities 

Diffuser performance predictions obtained from the inlet diffuser module 
are generally in good agreement with measured diffuser test data. The com- 
parison between predicted and measured pressure losses shown in Table I 
for three current GE combustion systems serves to illustrate the module's 
performance predictability. Greater accuracy in diffuser pressure loss 
predictions could be achieved by a more detailed modeling of the loss mech- 
anisms present. Airflow testing of full-scale annular models of the diffuser 
selected design is still necessary to provide final verfification of the pre- 
dicted pressure losses and evaluate the effects of distorted inlet velocity 


profiles and engine bleed flow requirements which are difficult to accurately 
analyze with potential flow models. 


Table I. Predicted Versus Measured Diffuser Ap/p. 





Outer 

Passage 

Center 

Passage 

Inner 

Passage 

GE29 Single Passage 

Predicted 

|3.11 

1.82 

2.84 

(M 3 =* 0.43) 

Measured 

l3.p8 

1,82 

2.85 

GE23 Single Passage 

Predicted 

12.59 

0.89 

2.45 

(M 3 =• 0.34) 

i 

Measured 

2.61 

0.74 

2.10 

CF6-80 Single Passage 

Predicted 

2.03 

0.64 

1.22 

(M 3 ■ 0.28) 

Measured 

2.06 

0.29 

0.89 


5.1.3 Airflow Distribution Module 
Description 

j The airflow distribution module (COBRA), is a computer analysis program 
that determines the one-dimensional aerothermodynaraic characteristics of a 
jet engine combustor. It provides a quick and economical method for computing 
the pressure and airflow distributions for a variety of combustor configura- 
tions. A flowchart of the COBRA program is presented in Figure 12. 

COBRA contains a general building block format which permits application 
flexibility without regard to the specific configuration of the combustor. 
COBRA is structured to contain a set of selectable options for the programmer 
so 1 that, as knowledge of aero thermodynamic systems is obtained, different 
types of combustor stations, additional types of pressure loss equations, 
and new types of airflow selections can be inserted rapidly into the existing 
program files. 

COBRA was developed using the GEBOSS Data Management System to take 
advantage of two of its main features: file storage and file editing. By 

using GEBOSS, all the source codes, data tables, and construction information 
are contained on one file called a Program File. This provides the user with 
a means of editing any of these files and greatly reduces the work of trans- 
ferring and maintaining COBRA at other installations. 



Figure 12^ Flow Diagram of Distribution Module. 







Analytical Methodology 

The COBRA program uses one-dimensional compressible flow relations. The 
basic procedure for an analysis of a combustor configuration with the COBRA 
program is characterized by the following steps: 

• Obtain a flowpath layout of the combustion system from the flowpath 
module and establish the station network as shown in Figure 13. 

• Identify the station type for each station in the model network. 
Definition of each station type are provided below: 

- Type 1; Normally used in region of combustor upstream of swirl 
cup or fuel injection plane. Options in the station allow for 
airflow selection, pressure drop selection, and hole flow. 

- Type 2: Liner cooling hole flow station designated in a non- 

burning pathway. Options are hole flow, airflow selection, 
apd pressure drop. 

Type 3: Flame pathway cooling flow station. Options are air- 

flow selection and pressure drop. Fuel-air ratio and tempera- 
ture rise are computed automatically. 

- Type 4: Liner dilution hole flow station designated in a non- 

burning pathway. Options are hole flow, airflow selection, and 
pressure drop. 

- Type 5: Flame pathway dilution hole station - inserted one 

slice downstream of companion skirt dilution hole flow station. 
Options include airflow selection, momentum pressure loss due 
to dilution flow entering flame passage, and other pressure 
drops. Fuel-air ratio and temperature rise are computed auto- 
matically. 

Type 6: Liner cooling and dilution hole flow station designated 

in a nonburning pathway. Options include separate cooling flow 
and dilution flow hole selections, airflow selection, and pres- 
sure drop. 

Type 7: Flame pathway cooling and dilution - inserted one slice 

downstream of companion skirt station. Options include airflow 

selection, momentum pressure loss due to dilution flow entering 
flame passage and other pressure drops. Fuel-air ratio and 
temperature rise are computed automatically. 

• Determine the required matrix for the stations, splices, and pathways 
(see Figure 14). These station types are assigned appropriate loca- 
tions on the flowpath layout and are therefore arranged in a pattern 
(rows). The program proceeds through the lattice by computing from 
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top to bottom on the slices. For program growth purposes, during 
the initial setup o£ lattice, slices may be placed through areas 
of the combustors where cooling or dilution stations may be 
needed at a future date. 

• For each station, select the proper data array options in the pro 
gram to determine airflow, pressure loss, and airflow splits. 
Options available for each station type are shown in Table II. 


Table II. Options Available for Types of Stations. 
(Options Designation Indicated by Circled Number) 


Station Type 

1 

2 

3 

4 

5 

6 

7 

Option 

Code 


1 AWC 


1 

AWD 


1 

AWC 




® 

W 

2 W 

1 W 

2 

W 

1 w 

2 

W 

1 

W 


@ 

TTO 

3 TTO 

2 TTO 

3 

TTO 

2 TTO 

3 

TTO 

2 

TTO 


® 

AP 

4 AP 

3 AP 

4 

AP 


4 

AP 

3 

AP 








3 APM 

5 

AWD 

4 

APM 









6 

AP* 




® 

AWFm 

5 AWFji 

4 AWFjj 

5 

AWFm 

4 AWFm 

7 

AWFm 

5 

AWFm 


® 

Rchw 

6 Rchw 

5 RcHW 

6 

Rchw 

5 RcHW 

8 

f^CHW 

6 

^^CHW 


© 

APfn 

7 APfn 

6 APpN 

7 

4Pfn 

6 APpN 

9 

APpN 

7 

APpN 


® 

APsp 

8 APgp 

7 APgp 

8 

APsp 

7 APgp 

10 

APsp 

8 

APsp 

Definitions: 











W 

- Weight flow calculation 

for station. 





WC & WD ^ Calculation of cooling or 

dilution hole 

flow. 



TTO 

" Total temperature option. 







P 

? Pressure 

drop calculation. 







PM 

■* Momentum 

loss calculation. 







WFm 

- Flow through hole, i.e. 

fuel nozzle hole in a cowl. 

i ^CHW 

- Recalculate hole 

flow (cooling or swirl 

cup 

, etc.) 



when J2 

or J6 station is a 

cooling path 

station adjacent 

i 

to 

initial flame 

path station. 

Used only on station 

1 

located 

on flame pathway. 







i PPN 

- Passage 

pressure 

drop due 

to a fuel nozzle 

stem or 



mounting 

pin. 









Psp 

- Passage 

pressure 

drop due 

to an 

airflow 

strut in 

passage. 


both profile drag and interference drag 

loss in calculated. 
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the airflovT of any other station or stations. In this way, flow splits can 
>e designated where desired and flow from two or more stations can be con- 
solidated. 

When liner cooling or dilution holes occur, the program determines the 
flow coefficient from an input curve of coefficient versus pressure parameter. 
Since there is room for a number of these curves, a liner hiving a number of 
different types of openings may be analyzed. These curves may be changed or 
replaced according to the user's discretion. 

Pressure loss considerations featured in the COBRA program include 
internal (one-dimensional) mixing losses due to momentum exchange, parasitic 
losses, and losses due to heat addition. These losses are represented by 
subroutines contained within the main program, and are automatically cal- 
culated depending on the type of station and option involved. In the case 
of parasitic losses, drag and area coefficients for the particular obstruc- 
tion must be included in the input data. In addition to these losses, a 
number of other pressure loss options are available through option input. 

These include: 

• Diffuser efficiency 

APt « (I-Hd) (Ps2 (IDEAL)-Ps^) (1) 

• Frictional losses 

Ap,j, = 4f (l/D)Q (2) 

• Coefficient of pressure (Cp). This option is usually used in a 
model when specific test data is available. 

• Percent of J8 total pressure. This option if primarily used to 
define pressure losses from station to station known from measured 
test data, but the Specific mechanism is not known. 

• Turning losses 

Ap^ = Q2(TURN) Qi (3) 

• Pressure loss due to ram across a carburetor cup 

^Tl|~ ^S2 ~ (Wj/Ae)^ (4) 

• Pressure loss due to pressure diferential across a swirl cup 
^Pg = R/2g (W^/Ae)^ 

PT2 * PS2 <Tt2/Ts2)^/y-1 

• Pressure losses due to screen or profiler placed across the flowpath. 


Percent of any station. This option is primarily used to specify 
a level of pressure at a given station. 


In performing the calculations, the program proceeds through the matrix 
and at each station performs the required calculations. A typical working 
layouts is shovm in Figure 15. The J1 through J8 designations indicate stations 
adjacent to the working station and represent an important key to the flexi- 
bility provided by COBRA (Figure 16). As an option, J1-J8 may be constructed 
automatically by the program, or this selection may be preempted by the user 
and arbitrary stations designated. 

COBRA uses a special subroutine to achieve the flow balance. This sub- 
routine is a simultaneous, multiparameter iteration technique. The balance 
package is reduced to a set of independent variables and error terms. The 
subroutine increments each of the independent variables and calculates the 
value for each of the error terras. The changes of the errors and independent 
variables are used to form a matrix of the partial derivatives. • This matrix 
is then used to calculate the variables at the next station. 


Prediction Capability 

The COBRA analysis program has been used extensively in the development of 
combustion systems for the CF6 engine family, the CFM56 engine family, and 
recently in the development of the NASA/GE Energy Efficient Engine (E^) com- 
bustor. Calculated flow parameters from models of these combustors have pro- 
vided accurate predictions which have been very useful in these development 
efforts. A comparison between the predicted airflow distribution from a 
COBRA model of the E^ double- annular development combustor Mod I configuration, 
and the airflow distribution calculated from measured data is presented in 
Table III. 


5.1.4 Combustor Heat Transfer Module 
Description 

I The combustor heat transfer module is comprised of two unlinked computer 
programs; HTTBC and SSHT. Both programs consider only steady-state heat 
transfer. The preliminary calculations are based on one-dimensional pro- 
cedures and are performed in HTTBC. This program is used along with the 
airflow distribution module in an iterative manner until the bulk liner metal 
temperature requirements are satisfied. Upon completion of this task, two- 
dimensional calculations are made to provide both maximum and nominal 2-D 
temperature distributions. This provides the required information to estimate 
the 3-D circumferential effect which is then input to the stress and life 
analysis procedures. These 2-D calculations are performed in SSHT. A flow- 
chart of the heat transfer module is provided in Figure 17. 















Table III. Predicted Airflow Distribution Versus Calculated Results 

from Measured Data, Double Annular Development Combustor, 
Mod-I Configuration, Flows Expressed in Terms of %w 
Combustor. 



Cobra Prediction 

From Test Data 

Outer Dome 

24.8 

25.4 

Inner Dome 

36.8 

38.3 

Centerbody 

CooHng 

5.8 

5.7 

Dilution 

5.9 

5.9 

Outer Liner 

Cooling 

6.7 

6.4 

Dilution 

3.4 

3.4 

Inner Liner 

Coo 1 ing 

7.2 

6.9 

Dilut ion 

9.4 

8.0 

Overall Total 

Pressure Drop 

5.5% 

6.0% 
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Figure 17. Flow Diagram of Liner Heat Transfer Module. 
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For transient calculations, the THTD program is used to calculate the 
time varying temperature distributions. For these transient calculations, 
heat transfer input parameters must be calculated for selected times in the 
transient using available test data or cycle deck calcultions. These values 
are then used as input along with the geometric description of the problem. 

The THTD program can perform three-dimensional analyses accounting for 
geometry, flow, and temperature features in the circumferential direction. 

The model preparation and calculation of boundary conditions is significantly 
greater than for two-dimensional analyses; hence, this approach is seldom 
used’ At General Electric it has been established that the three-dimensional 
effects of the details around coolant slot metering holes on the liner panels 
can usually be approximated satisfactorily with proper modifications to two- 
dimensional calculation. It was also determined that cyclic circumferential 
temperature patterns for thin liners involving one cycle per fuel nozzle, 
coiuld be adequately calculated with a series of two-dimensional calculations 
tiiat neglect circumferential conduction. However, thicker liners of the 
machined ring cooling slot type may have larger circumferential conduction 
effects which could require three-dimensionl analysis in support of similar 
stress and life calculations. 


Analytical Methodology 

I 

The initial total coolant flow and distribution is estimated from the 
flow required to cool similar or related designs for similar cycle conditions. 
The parameter used is the coolant flow per unit of cooled surface area per 
atmosphere of combustor pressure. General Electric has established values of 
this parameter as part of the combustion system design practice based on years 
of combustor development. This guide along with the total target coolant 
fiow is used to determine the initial coolant flow distribution. After the 
iiiitial flow distribution has been selected, the individual panel flows are 
examined to ensure that the levels are not less than a lower limit which has 
been established to maintain a protective film over the entire panel length. 
The flow distribution and combustor information are input to the Flow Distri- 
bution Module (cobra) which calculates the pressure distribution, the local 
gas velocity, fuel-air ratio, and gas temperature within the combustor. The 
velocity, fuel-air ratio, and gas temperature values are used in generating 
the required heat transfer input values. 

Both the one-dimensional and two-dimensional heat transfer analysis 
programs solve for a balance of the heat transfer quantities considered. 

These include flame radiation in hot side convection, conduction, backside 
(cold) convection, and radiation out. These heat transfer quantities are 
illustrated in Figure 18. In the case of the one-dimensional program, only 
conduction through the material thickness is considered. The results provide 
the calculated hot side and cold side metal surface temperatures. In the case 
of the 2-D program, solution to the heat balance is performed at each node of 
the entire node network defining the liner segment geometry. Since conduction 
in the axial sense is also considered, the calculation provides the axial 
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temperature distributions as well as temperature distributions through the 
thickness of the material. 

In each of these steady-state programs, correlations for heat transfer 
coefficients and film effectiveness are incorporated within the programs. 
However, it is necessary to use the existing General Electric data base as 
a means of determining which correlations to select as most appropriate to 
represent the design. The convective heat transfer coefficients are calculated 
from correlations developed from open literature data, and from wind tunnel 
test results for specific geometries. In the case of gas-side heat transfer, 
two correlations are available. One correlation is based on a wall jet flow 
directly downstream of the film slot and the second is based on a main gas 
stream correlation where the main stream velocity is relatively high or near 
the aft end of each panel where the wall jet velocity has decreased. Values 
are calculated from both correlations and the larger value is used in the 
analysis. The backside heat transfer coefficients are calculated from a 
cprrelation based on the flow of coolant between the liner and the casing or 
in the case of advanced liners the values are calculated from impingement 
flow correlations developed from model tests. The overhang is cooled by the 
impingement of coolant flow and the coefficient is calculated from a corre- 
lation based on model tests of the slot configuration, 

j film effectiveness data are based on wind tunnel test data of speci- 

fic film slot designs and on film temperature data derived from combustor 
tests. In areas where the film is interrupted, such as downstream of dilu- 
tion holes, the effective film must be back calculated from liner tempera- 
ture measurements. Normally the film effectiveness is calculated for a 
single panel ihd does not account for any film benefit from upstream panels. 

Hot streaks are simulated in the calculation by assuming circumferential 
profiles in velocity and temperature. 


The flame radiation levels are estimated from measured liner temperature 
data obtained on similar designs previously operated at similar conditions. 
This back calculation is done assuming a film effectiveness level and the 
local gas temperature calculated by the Airflow Distribution Module. The 
radiation level and film effectiveness combination is then used to analyze the 
design configuration. 

The basic procedure for a heat transfer analysis using the 2-D calcula- 
tljon approach is as follows. The first step is to lay out the node network 
or^ a 10 times enlargement of a cross section of the liner and digitize the 
cdordinates. The resulting x, y coordinate file is input to two time-sharing 
programs . One program (SSHTNHGN) prepares a node hookup file and a second 
prjogram (NPLTSSHT) prepares the x, y coordinate file in an SSHT main program 
fo^rmat. Other options are available in the NPLTSSHT program. One option pre- 
pares a plot of coordinate data which serves to check for possible errors in 
the model. Another option can scale the length of the model so that an exis- 
ting model can be used for designs where only the panel lengths have been 
changed. Next, the cycle and combustor geometry data time-sharing files are 
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^ fM geometry files, the node hookup files, and x, y 

coordinate file are merged to prepare a complete input file for the main 
program. The output from the main program is in the form of a hard copy of 

tSre"”^Th^°tiJf *=he calculated centroid temperl- 

ture. The time-sharing file is input to the NPLTSSHT program and prepares 
the temperature distribution plot. a ana prepares 

Predictive Capability 

j The liner heat transfer module has been applied to the design and analysis 
pf a wide variety of combustor designs for a number of years. T^se designs 

^ ^ Standard convection/ film cooled designs, such as the FlOl combus- 
tor, and more advanced impingement/ film cooled designs. The cooling slot 

from interrupted 360» slots, continuous 360* sJots s^ch 
as^the FlOl combustor and segmented film liner designs such as the e 3 combus- 

^ comparison between predicted axial temperature profiles and some 
measured metal temperature data for a FlOl liner panel is shown in Figure 19 
As observed from these comparisons, liner metal temperature predictions gener- 
ated from the heat transfer module generally reflect the measured levels? 

5.1.5 3-D Internal Flow Module 

the needed improvement in combustion system design and 
analysis, the combustor model currently in use at General Electric is beine 

to Include det.iUd 3-D interLl flow celculations! “Ser 

'“"prising the overall combuaeot model have sufficient sophistication 

eTfotr:ra -iot 

■ i 

a *3 flow module features three principle parts integrated within 

° program. The dilution region module is concerned with 

cooling air to guide exit gas temperature 
pattern development, including the hot streaks along the liners. The fuel 
insertion module is concerned with the initial specification of fuel droplet 
sxtze distribution, and velocities of fuel spatial distribution plus the 

vaporiS T^e^n b the flow field as the drops penetrate and 

kine^i«! <=‘'“bustion zone module handles heat release and chemical 

^ The initial purpose of the internal flow module will be to provide 

Drbvidrd^rM'‘‘'r5’”® of combustor exit gas temperature pattern tLnds, and 
provide detailed liner hot-side aerothermal properties to be used in con- 
junction with the heat transfer module for improved liner metal temperature 

plus iguition aSd°biomurchMM?eris?ic‘s?“ “'““Ohs. 
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Temperature 





The internal now module at General Electric features the INTFLOW computer 
program. It contains two alternate computer programs within its framework, 
one prepared by The Northern Research and Engineering Corporation, Reference 2 
and the second prepared by the Garrett Turbine Company, Reference 1. 

^ Due to the high degree of complexity in the program structure of these 
internal flow programs, discussions into their analytical methodology, and 
input/output parameters will be left for the reader to obtain from the refer- 
ences indicated. The analytical predictability of the internal flow module is 
to be assessed as the major portion of the Phase I Program effort. 

Some improvements were developed and incorporated into the 3-D internal 
flow model. This work was performed as part of an ongoing IR&D program at 
General Electric. The improvements include; 

• Direct individual node treatment of the swirl cup input 

• Grid refinement capability 

• Stair-step boundary treatment at the front and aft ends of the 
problem 

• Dilution injection on stair-step walls with injection angles 

• Simplified dilution injection input. 

In addition, an array of plotting programs were developed as part of the 
same IR&D program to provide visual interpretation of the output generated. 
Studies have been underway, however, on improved approaches for the basic 
calculation framework. As these new basic concepts are programmed. General 
Electric plans to replace the two programs. References 1 and 2, with an 
improved basic code. It will have features that have been identified in 
earlier studies at General Electric, from the efforts reported herein, and 
from the rapidly progressing technology expected to become available in the 
near future. 


2-D Turbulent Reacting Flow Codes 

To examine the details of the combusting flow, either experimentally or 
computationally, a geometrical simplification to two dimensions is required. 
The complex turbulence and turbulence-chemistry interactions of the fully 
3 d case along with the numerical difficulties are still present so that the 
2-D studies are very relevant. Furthermore, because fewer data points need 
to be taken and because the computations can be made on finer grids for the 
same or less storage and execution time than in 3-D, many more theory-data 
comparisons and of improved quality are available. 

The 2-p research code GETREF had already been developed at General 
Electric prior to this study. As discussed in Section 6.0, it contains the 


single~scalar assumed pdf/equilibrium chemistry/K-c turbulence model. Favre- 
averaging of the governing equation is employed exclusively. A parabolic 
fprmulation is used to simulate shear flows, for example, a burning jet in 
cpflowing air, and uses a forward-marching solution algorithm. The elliptic 
version involves a standard relaxation technique: the primary numerical 

inaccuracy is that caused by convection term discretization. These points are 
discussed in much greater detail in Section 6.0. In these codes, isothermal 
flows are simulated by suppressing the combustion aspects; this is facilitated 
by their modular construction. 


5.1.6 Usefulness of Output Parameters 

The relative importance of the various modules comprising the overall 
combustor aerothermal anaysis model are listed in the order of usefulness in 
General Electric Design procedures below: 


Airflow Distribution Module 

• Predicted external and internal velocities and pressures 

• Predicted pressure drops 

• Predicted flow distribution or flow areas. 

This output is considered most useful because it is used on a daily basis in 
the design and development of combustor?. The impact of design modifications 
on airflow distribution and pressure drop can easily be analyzed. 


Heat Transfer Module 

The metal temperature distribution is used by stress analysis programs for 
estimating life and is used for evaluating the effect of configuration changes 
(primarily cooling changes) and postulated changes. 


Internal Flow Module 


Although not currently in design use, INTFLOW potentially has the ability 
to predict velocity flow fields near the liner and perhaps hot streaks with 
improved accuracy that will greatly improve the accuracy of the heat transfer 
module output. It is also expected to provide useful guidance as a pertur- 
bation calculation for dilution hole trimming of pattern factors. 


Inlet Diffuser Module 

The design technique minimizes development time and provides total 
pressure levels in each of the passages around the combustor for input into 
the Airflow Distribution Module. 
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Preliminary Design Module 


Combustor sizing parameters provides envelope that generally remains 
unchanged throughout development work. Similar information can, however, come 
from existing correlations from many sources other than this particular com- 
puter program. 

Light-off, Blowout, Combustion Efficiency and Emissions 

These parameters are currently estimated by empirical correlations with- 
out using the calculation modules. 


5.2 S ELECTION OF BENCHMARK QUALITY DATA BASE 

One of the specific objectives of the Phase I effort is to identify defi-’ 
ciencies in the aerothermal models currently used in the design and analysis 
of gas turbine combustors. Towards this purpose, a data base was assembled 
and studied for selection of benchmark quality test cases used in evaluating 
each model in the overall combustor aerothermal model. This data base includes 
data from the available literature, existing General Electric combustor/ engine 
test data, and additional new data obtained from the experimental testing 
effort conducted as part of the overall Phase I modeling program. 


5.2,1 Available Benchmark Experiments from the Literature 

I The flows which occur in the combustors of gas turbines involve varia- 
tions of mean properties in three directions with complex turbulence charac- 
teristics, combustion which is largely of the diffusion-flame type but may 
involve some burning of premixed fuel and air, and the evaporation and combus- 
tion of liquid droplets. Calculation methods, based on the solution of con- 
servation equations in differential form, require the adoption of models for 
these physical phenomens as well as schemes to transform the differential 
equations into finite-difference equations. Thus, even the exact solution to 
t^e finite difference equations which are actually solved by approximate 
iterative procedures) would not exactly represent the flow. The validity of 
the models can, in principle, be evaluated directly or indirectly. Thus, for 
example, the components of an effective-viscosity assumption. 



can all be measured and the value of pgff can be determined explicitly. In 
general, however, the direct appraisal of assumptions is not possible and the 
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merits of all, or most, of the assumptions in a calculation method must be 
appraised by comparison of calculated and measured results. 

While only a few experiments from the literature were selected for com- 
parison with data in this program, a somewhat broader group of experiments 
were considered. The discussion that follows while not an exhaustive litera- 
ture survey does indicate many of the available experiments and very briefly 
presents some of the significant findings from comparisons with models that 
have been previously emphasized in the literature. 


Isothermal Flow 

In the past, many experimental and computational investigations have been 
carried out with two-dimensional flows, that is, those which can be described 
by two independent spatial variables. The confined-flow measurements of, for 
example, Pai, et al., 1975H; Owen, 197612; a^d Habib and Whitelaw, 198013, 
come into this category and have allowed calculation methods, such as that of 
Mongia, 19791, to be tested against sudden-expansion configurations and in the 
presence of swirl. It is clear from the various related computational inves- 
tigations that the careful location of a large number of grid nodes is required, 
even in two-dimensional flow, to avoid serious numerical errors. In addition, 
the advantages and disadvantages of different turbulence-model assumptions have 
been demonstrated and it can be concluded, for example, that a transport model 
is necessary to represent recirculating flows and that high degrees of swirl 
are not well represented by isotropic eddy-viscosity models such as the two- 
equation forms proposed by Jones and Launder, 19721^, or Saffraan, 197415. At 
the same time, the use of numerical grids based on, for example, 400 nodes in 
a two-dimensional solution domain, will usually give rise to cell Reynolds 
numbers greater than |2| generally leading to first order upwind differencing 
and numerical diffusion. The introduction of numerical diffusion often nulli- 
fies any added benefit which may be gained from the use of Reynolds stress 
models. Since the problem is even more acute in three-dimensional flows, it 
is questionable whether higher order closure in turbulence models (such as 
Reynolds transport) can be justified for use in real combustors. Furthermore, 
turbulence exerts less control in the jet-in-crossflow arrangement of combus- 
tbrs, where pressure forces are important, than in the sudden— expansion, 
swirling flows mentioned above. On the other hand, in the future it may prove 
necessary to further investigate the relative merits of turbulence models in 
two-dimensional, turbulence-controlled flows; however, with practical calcula- 
tion methods required in the near future, the present investigation will con- 
centrate on flows which possess features more directly relevant to gas tur- 
bines . 

The complex geometry of a gas-turbine combustor, even in the absence of 
combustion, involves the interaction of primary jets, dilution jets, the flow 
from fueling devices and film-cooling jets in the presence of an annular or 
cannular confinement. A simplified axisymmetric arrangement was investigated 
by Green and Whitelaw, 19805, This is one of the experiments selected for 
comparison with calculations in this program. 


The experiment features an axisymmetric form of the crucial dilution jet 
problem. These data allow examination of the prediction with emphasis on 
numerical accuracy and the resolution of pressure gradients and without the 
problems of chemical reaction. 

The data of Vu and Gouldin, 1980^ for a double-annular confined flow with 
both jets CO- or counterswirled are also an excellent data set for comparison. 
Detailed velocities and turbulence quantities were obtained over a range of 
swirl numbers. Inlet quantities were clearly measured, thereby providing a 
benchmark quality data set. This particular data was used in the present 
efforts to test a modified form of the two— equation model (the algebraic 
stress model) which directly accounts for the anisotropies of eddy viscosity 
caused by the streamline curvature. ' 


The round jet in crossflow represents another simple flow used for model 
testing purposes. It possesses the basic features of one important aspect of 
combustor flows and permits determination of the influence of node density and 
distribution particularly in the near jet-exit region. Measurements exist for 
a range of velocity ratios. Table IV provides a list of the available data. 

It should be noted that previous computational investigations of the round jet 

Patankar, et al. 197316; jones and McGuirk, 
1979 , Crabb, 1979,1981 » , and White, 1980^0, and provide guidance. Never- 

theless, this well-measured flow provides a foundation upon which any calcula- 
tion development can be based. One particular point of note from previous 
calculations is the conclusion of White that a greater distribution of nodes 
than his 25 x 20 x 15 is required to achieve grid independence. A careful 
analysis of the numerical truncation error, Claus, 198321, indicates that at 
least 75,000 nodes are needed. 


The next level of complexity is the row of jets in crossflow and opposed 
rows of jets in crossflow. In the former case, measurements have been pro- 
vided by Kamotani and Greber, 197222; Holdeman and Walker, 1977^; Crabb, 

1979 , and Khan and Whitelaw, 198023, accuracy of measurements made 

mainly by pitot tube is poorer than that of the single jet data but is more 
than adequate to assess the calculated mixing of the jets as a function of 
pitch to diameter ratio. It should be noted that the finite and known width 
of the solution domain renders this flow, in this one sense, easier to repre- 
sent by a computational mesh than the single jet in crossflow. Again, related 
calculations have been reported, Khan, et al., 198224, and provide useful 
guidance but in no way obviate the need for further calculations. 

The data of Holdeman and Walker 19779, and also reported in Walker and 
Kors 1973 , were selected for comparisons in this work partly because these 
data have been previously used at General Electric in guiding combustor design/ 
development efforts. 

The flow downstream of opposed jets has been measured by Kamotani and 
Greber, 1974'^^ and Atkinson, et al., 198226 and, to date, no related cal- 
culations have been reported. The complex structures reported by Atkinson, et 
al.-^o, and created by pressure rather than turbulence forces could provide a 
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Table IV. Round Jet in Crossflow 
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Table IV. Round Jet In Crossflow (Concluded) 
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where P(f) represents the density-weighted probability density distribution of 
scalar property f. In the Raman scattering data, Section 7.1, both types of 
average were derived from the instantaneous data taken. 

A very large number of gaseous-flame measurements is available in the 
literature ranging from open diffusion flames to confined premixed flames and 
to real combustors. Again, a careful choice was made to concentrate on the 
confined nonpremixed configuration. For completeness, however. Table V oro- 

I f “E avilf U blu(£-body premlied aid Sr dltr 

which may also be used for testing models. 

In calculating combustion flows, one must realize that the several models 
combuitif "n numerical solution techniques, turbulence, heat transfer, and 

uncertainties in proportions and with magnitudes that 
particular flow being investigated. Problems with 
nijmerical diffusion due to coarse grids have been discussed, and two-equation 
turbulence model weaknesses in, for example, recirculation zones have been 

II dd.fr'”"® • ■ continued improvements in these areas 

Ilr-A H comments can be made regarding 

kinetics in the combustion model, the need for handling complex 

“ E"" of complex fUm-cooled 

I The calculations of Jones, 197928 have shown that the hydrogen-diffusion 
fl^ame of Kent and BUger, 197329 represented with fast chemistry 

assumptions and that different assumed forms of the scslar probability^density 
function can lead to different results. General Electric has just completed a 
Department of Energy-funded study, Lapp, et al., 1983^ on turbulent jet dif- 

mIaI°LrfT^^h on mean and fluctuating velocity and Raman data on 

mean and fluctuating composition and temperature are thus available. This 

tillTe turbulence/chemistry model in a computationally 

pi:!e:te]°rb sec^”: 7af‘ “"">“lso„. and fha caauUa ana 

man dj-ffusion-flame flows include those of Owen, 197612; 

ffIpH El liman, et al., 197831. In both cases, the flow is con- 

thl cl,rnf^v''^‘'’30"3l“r^ ® hydrocarbon fuel in a nonpremixed manner. In 
the case of Elliman,3 ,31 example, nonreacting and reacting flow fields 

thrl annular baffle with methane flowing 

k f'?!' ‘^^^“^ter central hole and air flowing with 12.5 mm space 
iTrnrT. confining 150 mm diameter pipe. Measurements of 
velocity and temperature were obtained with a seven-hole spherical impact 

Sled^itl II thermocouple. The local equivalence ratio was meas- 

ured with an uncooled i 4mpling probe. The measured values correspond to 
density-weighted averages and the uncertainties in velocity (outside the region 

III" temperature, and equivalence ratio are not greater 
phan 10/^, 15^, and 20%, respectively. 

The flow investigated by Toral and Whitelaw, 198232 £„ a sector rig 

complexities associated with a semicircular combustor head ’and 
film cooling jets which drive the primary vortex in addition to cooling the 
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Table V. Premixed Combustion (Concluded). 
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wall in critical areas. The calculation of this flow requires preliminary 
steps toward a general geometry package and could provide a very helpful 
building brick in the development of a method which is generally applicable to 
practical combustor arrangements. This is particularly so since the data 
include isothermal measurements and distributions of a passive scalar with and 
without combustion. Thus, the numerical and turbulence model can be appraised 
itii relation to the isothermal flow velocity results; the heat transfer model 
can be evaluated in relation to the isothermal flow scalar distributions; and 
the combustion model can be examined in relation to measurements of velocity, 
temperature, and the concentrations of unburned hydrocarbon, carbon monoxide, 
carbon dioxide, and oxygen. It is useful to note that the authors deduce from 
the measurements that equilibrium or fast chemistry assumptions appear to be 
justified except in rich regions of their flame. 

The axisymmetric bluff body diffusion flame data by Lightman, et al. 

1980^ is another source of detailed velocity data acquired with an LDV system 
for both nonreacting and reacting flow situations over a broad range of both 
fuel and air flows. This experiment was modeled in the present study and will 
be discussed in detail (Section 7.1). 

Less detailed results, in combustor geometries which are even closer to 
practice, have been reported for example by Tuttle, et al., 1973^3; Vranos and 
Tabak, 19763^; and Noyce, et al., 1981^5. 


Droplet Combustion 

Models for droplet combustion are usually based on Crowe, 197436 and 
later publications, for example that of Crowe, et al.. 197737, They can allow 
for the sCffaration of the droplet and gaseous velocities by a Lagrangian track- 
ing procedure or can assume no slip between the two phases. Which is more 
applicable depends on the problem being considered. With the levels of pre- 
heat temperature present in many gas turbines and the atomization of the fuel 
td droplets with Sauter mean diameters of less 50 Pm, it is unlikely that the 
cdnsideration of slip is required, except possibly at idle conditions. In 
addition, the assumptions necessary in the numerical, turbulence, heat trans- 
fer, and combustion models suggest that the droplet model should be as simple 
as possible. Appropriate assumptions for evaporation and group effects can 
readily be incorporated as required. 

The application of the procedures of Crowe have been reported by, for 
example, El Banhawy, et al., 1978,38 El Banhawy and Whitelaw, 1980^^ and for 
comparison purposes have made use of the experimental data of Tuttle, et al., 
1973^3; These data correspond to a hollow spray of monodisperse kerosene drop- 
lets of around 50 pm with no preheat. It appears from these and similar 
results that the droplet tracking procedure will only be necessary for the 
idle condition but further evaluations could be carried out in relation to the 
above data and to that provided by studies such as those of Botros. et al. 
1980^0; Driscoll and Pel accio, 1980^1; and Nicholls, et al., 1980^^. 


Experiments Used to Assess Model 


.u summarizes the experiments selected along with the type of data 

gathered. Again the two-dimensional elliptic code GETREF and a parabolic 
version of it were used for the axisymmetric experiments to assess the mathe- 
matical/physical model. For the Walker and Kors, 19738 ^ata, calculations 
were accomplished with the 3-D INTFLOW code which provided an overall assess- 
ment of the code though with a much less detailed grid and without the error 
analysis performed for the axisymmetric assessments. 


General Electric Engine and Component Data Base 

VI Electric, large quantities of measured test data are avail- 

able on all production and development combustion systems. These data in- 
clude diffuser performance measurements, combustor airflow distribution/hole 
flow coefficient data, combustor liner temperature data, and combustor exit 
gas temperature pattern data. In addition, measured data are available for 
combustor stability and exhaust emission performance. This extensive data 
base has been instrumental in the establishment and refinement of combus- 
tion systems design guidelines used at General Electric. 

j performance data are routinely obtained from full-scale annular 

model testing of each new diffuser design. This testing is done as a verifi- 
cation of the design performance predictions obtained from the inlet diffuser 
module. Detailed pressure measurements are made on numerous positions along 
the diffuser model flowpath. Such pressure data are obtained over a wide 
range of operating conditions providing sufficient information to generate 
performance design curves. These curves become very useful to the designer 
in predicting changes in diffuser performance at off-design operating condi- 
tions. A schematic of a typical model test rig, and performance curves 
are presented for the NASA/GE Energy Efficient Engine diffuser design in 
Figures 20 through 23. Because of the finite amount of instrumentation used 
ip these tests, specific details concerning circumferential nonuniformity, 
boundary layer effects, and wakes downstream of support struts are not well 

the vast amount of design experience accumulated has 
to the establishment of design guidelines sufficient to provide diffuser 
designs that demonstrate excellent performance. 

I A large volume of data concerning flow coefficients and bulk injection 
angles for different shaped holes has been accumulated at General Electric 
oyer many years. These data were obtained from detailed model wind tunnel 
which a large array of hole shapes were tested to define the flow 
coefficient and bulk injection angle. Data of this nature are available for 

General Electric combustor designs. Much of 
thas dilution hole data have been empirically correlated and introduced into 
current design methods such as the airflow distribution module. The impor- 
tance of having this data base readily available is illustrated in Figure 24 
by the large variation in hole flow coefficient that occurs as a result of a 
relatively minor change in the hole shape. These data also take into account 


Table VI. Experiments From the Literature Which Were Modeled. 


Isothermal Combustor-Like Flows 

1. A. Green and J.H. Whitelaw, Measurements and Calculations of 
Isothermal Flow in Axisymmetric Models of Combustor Geometries, 
Journal Mech. Eng. Sci., Vol. 22, No. 3, 1980. (LDV Measurements of 
Mean Velocity and Turbulence in Slit-Jet/Plenum Chamber Geometry). 

2. B.T. Vu and F. Gouldin, Flow Measurements in a Model Swirl Combustor, 
AIAA Paper No. 80-007, 1980. (Axisymmetric double-annular flow geom- 
etry where both flows were independently swirled. Mean velocity and 
turbulence measurements.) 

Combustion Flows 

1. M. Lapp, M.C. Drake, C.M. Penney, R.W. Pitz and S.M. Correa, "Turbu- 
lent Combustion Experiments and Modeling." Final report DoE con- 
tract DE-AC04-78 ET 13146. (Turbulent CO/H2/N2 jet diffusion flame 
in coflowing air. LDV and Pulsed Raman data for mean and fluctuating 
components of velocity, concentrations and temperature. Pdf's 
directly measured.) 

2. A.J. Lightman, et al, "Velocity Measurements in a Bluff-Body Diffu- 
sion Flame." AIAA Paper No. 80-1544. (Non-premixed axisymmetric 
bluff“body stabilized propane flame, LDV measurements for mean and 
fluctuating axial velocities. Thermocouple laeasurements for tempera- 
ture. 

3-D Dilution Jets 

1. R.E. Walker and D.L. Kors, "Multiple Jet Study Final Report". NASA 
CR-121217. Three-dimensional dilution jet experiments. Thermocouple 
measurements of jet penetration and mixing. 
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the effect of back side velocity and pressure drop on the flow coefficient and 
th^ bulk injection angle. As a verification, the flow areas of the individual 
flbw features of combustors are periodically measured in an atmospheric cali** 
brated flow stand. Despite the extent of the data base, there is a lack of 
data concerning turbulence levels in and around the penetrating jet, metastable 
jet characteristics, and details within the jet such as velocity and injection 
angle variation across the jet as a result of the nonuniform pressure fields 
surrounding the jet. These are not so important for the gross analysis con* 
ducted using the airflow distribution module. However, their importance is 
evident when attempting accurate detailed internal flow calculations. 

Measured combustor liner temperature data have been accumulated at 
General Electric for a variety of designs, both production and devleopment. 
These include conventional convection-cooled designs, impingement-cooled 
deiigns, and designs featuring convection and impingement cooling. This 
exfcei.s U'S data base includes the effects of fuel type, inlet pressure, inlet 
teiiiperature, and fuel/air ratio. A typical example of these data is provided 
in , Figure 25. The majority of these data have been obtained from component 
testing in which test conditions are set to simulated actual engine operations. 
Component testing provides the advantages of using large quantities of thermo- 
couple instrumentation, plus the ability to investigate a wide range of opera- 
ting conditions, both on cycle and off cycle. Additional temperature data are 
obtained from actual engine tests. The increased complexities of engine test- 
ing requires the use of significantly leas instrumentation, plus the inability 
to set operating conditions very far off the engine operating line. The advan- 
tages of obtaining liner temperature data from engine testing is that actual 
engine cycle conditions are set at high power settings, plus one has the abil- 
ity of oi)taining transient data during throttle movements. Figure 26 provides 
an I example of the transient liner temperature data base, these data still 
represent finite discrete measurements. As such, this lacks the ability to 
provide detailed information on axial and circumferential temperature varia- 
tion and detailed information concerning liner temperature patterns around 
dilution holes. At General Electric, it has been established that cyclic 
circumferential temperature patterns for thin liner shells involving one cycle 
per fuel nozzle can be adequately represented using this data base along with 
a series of two dimensional calculations. However, liner designs with thicker 
shells may have larger circumferential conduction effects which may require 
three dimensional treatment. 

In addition to liner temperature data, a substantial amount of cooling 
slbt film effectiveness data exists at General Electric. These data are based 
on wind tunnel testing of models of various cooling slot designs, and on film 
temperature data derived from combustor tests. These wind tunnel model tests 
cannot simulate combustion conditions such as swirling and recirculating flows, 
but are very useful in separating design features of various slots. These 
features include: 

• Feed hole spacing to diameter ratios. 


Feed hole orientation. 
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Figure 25. Sample of Steady State Liner Metal Temperature. 
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Figure 26. Sample of Transient Liner Metal Temperature Data. 




• Cooling slot overhang length. 

• Overhang thickness to slot height ratio. 

• Slot flow to mainstream flow ratio. 

• Slot plenum volume. 

Film effectiveness levels in actual combustor environments are derived 
from component or engine test data. These levels are derived either by back 
calculations from measured metal temperature data dr from direct measurement 
of cooling film temperatures. 

The General Electric Company has long recognized the importance of de~ 
veloping combustion systems which demonstrate ^aw levels of pattern factors 
and averaged profiles. The pattern factor paraifiteter is defined as: 

“ ^MAX 4 ~ ^3 

TavG 4 - T 3 (9 


where; 


TmAX 4 is the highest temperature anywhere in the exit plane, 
Tavg 4 the average temperature in the exit plane. 


is the compressor discharge temperature. 


The average profile is defined by: 
" ^LOCAL 

Tavg 


i' J 1 

i 1 1 -■ 


where ; 


“^LOCAL “ '^AVG 4 one radial location) - T^vq 4 
TaVG = 4 “ T3 


n h--: 


I The pattern factor defines the magnitude of the highest temperature level 
in the combustor exhaust stream, and is of concern to the structural integrity 
pf the stationary turbine hardware, particularly the turbine inlet nozzle 
Ijardware. The averaged profile is a measure of the circumferentially averaged 
radial temperature profile, and is of concern to the turbine rotating hardware. 
In support of these concerns, a substantial amount of exit gas temperature data 
have been obtained at General Electric. The majority of the data are acquired 
through detailed exit gas temperature measurements of full annular combustors. 






Some additional data are provided from exit gas temperature 
sector combustors. Sector data are cheaper to obtain, but can 
effects In all current as well as in past production programs, it has been 
t^I p«;tice to measure exit gas temperature patterns on production combustors 
on a^^sampling basis. This is done to ensure that profile and pattern facto 
levels a?e bling maintained in the production engines. A significant amount of 
additional exit gas temperature data are obtained from combustor designs in the 
development phase. This is done to ensure that the hardware ^ 

production sLisfies the design intent. So extensive is this data base, that 
Lpirical design information has been generated for all ^ 

which can be used to predict the impact on the profile and pattern factor from 
changes in airflow distribution, changes in dilution patterns, changes in fuel 
type and injection hardware, as well as the impact of 56 

deterioration. Table VII summarizes the results obtained 

combustors subjected to various hours of rigorous endurance ?' factor 

orkuction situation some variation in the measured profile and pattern factor 
Lh be encountered due to the mechanical tolerances allowed on the various 

combustor flow features and leakage. Sufficient data ^ An 

statisti-al experience range for the measured profile and pattern factor. An 
examptr;f this is shown in the following tabulation of data for a member of 
CFM56 combustors of the same design configuration. 


No. of 
Combustors 


Pattern 

Factor 


1 

3 

1 

1 

1 


1.24 

1.22 

1.21 

1.18 

1.17 


The exit gas temperature data base also serves to illustrate another 
ooint Typically there is a considerable amount of nonuniformity in the 
exifgas tempJraLre pattern around the exit annulus of a combustor. Figure 
27 shLs contours of the pattern factor parameter calculated from measured 
exit gas temperature data obtained from a typical FlOl combustor. 
observed from this figure that similarities exist from ®'^P ’ 

the variations in levels from cup to cup serve to illustrate the degree of 
nbnuniformity. This is of particular interest when considering three 
dimensional internal flow calculations. Current aerothermal ^ 

type typical lly feature cyclic boundary conditions. This essentially states 
thL the combustor section modeled is exactly identical to all other sections 
comprising the entire combustor. With current grid stations, models _ 
typically represent single cup sections of the entire combustor in order to 


71 




Table VII. FlOl Combustor Performance Stability. 


Combustor 

Pre-engine 

Pattern 

Factor 

Post-engine 

Pattern 

Factor 

Test 

Hours 


Comment 
ight Damage 

Pre-PFRT 

1.173 

1. 180 

lo 

! 

SI: 

PFRT 

1.215 

;.230 

1^2 

No 

Damage 

PFRT Development 

1.218 

i .255 

480 

Slight Damage 

PV 

1.220 

1.220 

390 

1 

Slight Damage 

Development 

i 

1.265 

1.263 

349 

Very Damaged 

PVj 

1.200 

1.210 

373 

Some Damage 

Development 

1.230 

1.200 

119 

Some Damage 

Development 

1.230 

1.260 

392 

Minor Damage 

Initial Product 

1.220 

1.200 

522 

Minor Damage 
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provide the grid detail necessary for accurate flow calculations. Thus, the 
current aerothermal model cannot address the circumferential nonuniformity 
issue. This represents a significant limitation in the application of the 
model for predictions of pattern factor. 

As previously mentioned, additional test data concerning combustor 
stability and gaseous and particulate emissions are also available at General 
Electric. Because they are not of significance in the scope of the Phase I 
effort, no further discussions about these data will be made. 

In selecting which sets of data to use as the benchmark test case for 
comparison with calculations from the General Electric aerothermal model in 
general use, consideration was given to a number of items. First, the com- 
bustion system represented by the data had to be of advanced, modern design. 
Additional considerations included the amount, quality, and availability of 
the test data, as well as the range of operating conditions represented in 
the data. Consideration was also given to the availability of wind tunnel 
model test data to document hole flow coefficients, and film effectiveness 
and heat transfer coefficient data. From these considerations, the data 
base for the GE/FlOl combustor was selected as the benchmark test case. 


5.3 EXPERIMENTAL TEST PROGRAM 

j As part of the Phase I Aerothermal Modeling Program, an experimental 

test program was conducted at General Electric. The purpose of this effort 
was to obtain additional definitive data with respect to the internal flow 
field of a combustor from which to assess the capabilities of the 3-D elliptic 
detailed internal flow model. A total of 18 configurations were tested using 
two available test rigs. Each test rig represented a 90“ sector of a full 
annulus combustor. 

The majority of the experim-;ntal testing was conducted in a rig repre- 
sehting a simple combustor design. This rig featured a flat uncooled dome 
inlet, and parallel cylindrical uncooled liner walls. The geometry of this 
ri^ was well within the geometrical constraints of the 3-D elliptic model. 

The dome section was interchangable to allow various dome configurations to 
be tested. In this testing effort four dome inlet configurations were tested 
in this rig. These configurations were: 

• Low pressure drop urtiform/parallel flow inlet using a wire mesh 

• High pressure drop uniform/parallel flow inlet using a perforated 
plate 

• Nonswirl air injector cups 

• Actual GE/FlOl swirl cup hardware. 

I 

The liner walls featured two rows of dilution holes patterned after the GE/FlOl 
combustor. The hole sizes could be varied or entirely blocked off depending 


on the desired configuration. Since this rig was totally uncooled, only two- 
temperature trace experimental testing was conducted using this rig. 

The parallel wall test rig was designed with three separate flow control 
systems to independently set conditions feeding the dome inlet, the outer 
dilution holes, and the inner dilution holes. This provided the capability for 
two-temperature trace testing. 


Data were obtained by using a rake with 21 equally (radial) spaced 
thermocouple elements mounted onto a traversing arm designed to allow axial as 
well as circumferential movement of the rake. During a typical test run, the 
rake was moved from -18“ from top center to +18* from top center at 1“ incre- 
ments at 4 axial planes downstream of the aftmost dilution injection point. 
This provided an array of 777 temperature measurements within the 36* section 
at each axial plane. In addition to the thermocouples, 11 equally (radial) 
spaced impact pressure probes were initially mounted onto the rake. These 
pressures were used to map out the inlet velocity profiles in the rig when 
using the uni form/ parallel flow dome inlet configurations. All data obtained 
was sent through a data logger system and loaded onto a cassett tape for 
posttest processing. 


A schematic of the parallel wall test rig is shown in Figure 28. An 
overall view of the test rig hardware and closeups views of various features 
of the rig are shown in the photographs in Figures 29 through 32. 

The other test rig used in the experimental test effort was a 90" sector 
of Ian actual GE/FlOl combustor. The FlOl combustor features a conical dome 
with 20 equally spaced counterrotating swirl cup assemblies. The liners are 
a machined ring type featuring 6 film cooling slots in the outer liner, and 5 
film cooling slots in the inner liner. The liners are contoured resulting in 
a Converging internal flowpath. Each liner has 2 rows of dilution holes. 
Corresponding rows on each liner feature opposed holes of equal size. The 
primary injection row features 40 holes equally spaced around the circumfer- 
ence, one hole directly in line with the cup center, and one hole between cup 
centers. The holes in line have a diameter of 0.84 cm (0.33 inch), while the 
holes between are larger with a diameter 1.04 cm (0.41 inch). The secondary 
injection row features 60 holes equally spaced around the circumference. The 
holes are offset 3* from the cup centers. All 60 holes have diameter of 0.91 
cm (0.36 inch). 


With this test rig, data were obtained by using a rake with 13 equally 
(radial) spaced thermocouple elements. The rake was mounted onto a traversing 
arm that permitted movement only in the circumferential sense. Thus tempera- 
ture measurements in this rig were made only at the discharge plane. During 
a typical test run the rake was moved from -36* from top center to +36* from 
top center in 1.5* increments. This provided an array of 637 temperature 
measurements within the 72* section of the discharge plane. All data obtained 
were sent through a data logger system and loaded onto cassett tapes for post- 
test processing. An illustration of the FlOl sector combustor test rig is 
shown in Figure 33. 
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Figure 28. Aerothermal Modeling Experimental Test Rig (Unfueled) 
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Of the 18 experimental configurations tested, 14 were conducted as two- 
temperature trace experiments in the parallel-wall test rig. The remaining 4 
experimental configurations were tested in the FlOl sector combustor test rig. 
Of these, two were conducted as two- temperature trace experiments ana two were 
conducted with fuel and heat release. This experimental test series started 
with the simple case of a single row of uniform jets penetrating into a 
uni form/ parallel crossflow, and progressed in complexity to cases in actual 
combustor hardware with fuel injection and heat release. This approach per- 
mitted investigating the effects of such phenomena as opposing rowi of dilu- 
tion jets the introduction of swirl at the dome inlet, actual flowpath con- 
tours and the introduction of fuel and heat release. A summary of the experi- 
mental configurations is provided in Table VIII. The test data obtained 
from this program provided an excellent source from which to make comparisons 
with calculations of the experimental configurations performed on the 3-D 
elliptic detailed internal flow model. 

All of the test data obtained are contained along with a detailed 
description of each experimental configuration and the test conditions in 
olume II. Also included in this volume are contour plots of the tempera- 
ure difference ratio for all two-temperature trace experimental test runs, and 
contour plots of the pattern factor parameter for all experimental test runs 
with fuel and heat release. The definitions of these nondimens ional tempera- 
ture parameters are provided below. ^ 

Temperature Difference Ratio: 


^Inlet ~ ^Local 
'^Inlet ■ Tjex 


( 11 ) 


Pattern Factor 


^Local ~ ^Inlet 

TavG - Tinltc 


( 12 ) 


Som» additional measurements were made in support of the overall experi- 
mental testing effort. Part of these additional measurements involved obtain- 
ing fuel spray patternation data. The other part involved making velocity 
measurements of the airflow field at the discharge of the FlOl swirl cup 
assembly. These data provided the boundary inputs necessary to model the dome 
inlet conditions for experimental test configurations featuring this hardware. 


test 


The fuel spray patternation data were obtained in the 
stand at General Electric. This test stand is shown 


radial patternation 
in Figure 34. In 
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Of the 18 experimental configurations tested, 14 were conducted as two- 

parallel-wall test rig. The remaining 4 
experimental configurations were tested in the FlOl sector combustor test rie. 
Of these, two were conducted as two-temperature trace experiments ana two were 
conducted with fuel and heat release. This experimental test series started 
with the simple case of a single row of uniform jets penetrating into a 
uni form/ parallel crossflow, and progressed in complexity to cases in actual 
combustor hardware with fuel injection and heat release. This approach per- 
mitted investigating the effects of such phenomena as opposing rows of dilu- 
tion jets, the introduction of swirl at the dome inlet, actual flowpath con- 
tours, and the introduction of fuel and heat release. A summary of the experi- 
mental configurations is provided in Table VIII. The test data obtained 
from this program provided an excellent source from which to make comparisons 
with calculations of the experimental configurations performed on the 3-D 
elliptic detailed internal flow model. 

All of the test data obtained are contained along with a detailed 
description of each experimental configuration and the test conditions in 
Volume II. Also included in this volume are contour plots of the tempera- 
ure difference ratio for all two-temperature trace experimental test runs, and 
contour plots of the pattern factor parameter for all experimental test runs 
with fuel and heat release. The definitions of these nondimens ional tempera- 
ture parameters are provided below. 

Temperature Difference Ratio; 


0 3 , ^Inlet ~ ^Local 

'^Inlet " Tjex 


( 11 ) 


Pattern Factor 


Pf 
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Soma additional measurements were made in support of the overall experi- 
mental testing effort. Part of these additional measurements involved obtain- 
ing fuel spray patternation data. The other part involved making velocity 
measurements of the airflow field at the discharge of the FlOl swirl cup 
assembly. These data provided the boundary inputs necessary to model the dome 
inlet conditions for experimental test configurations featuring this hardware. 

patternation data were obtained in the radial patternation 
test stand at General Electric. This test stand is shown in Figure 34. In 
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this test stand, the swirl cup - fuel injector assembly is mounted centrally 
above a multitude collector rake. The tubes in this rake are arranged along a 
radial arc and serve to collect fuel over a predetermined run time. The plate 
onto which the swirl cup is mounted can be rotated with respect to the rake 
allowing measurements of radial fuel distributions at any circumferential 
position. Test runs are conducted at atmospheric pressure and ambient temper- 
ature. The swirl cup and fuel injector pressure drops can be set to simulate 
the desired test rig or engine operating conditions. Radial patternations 
were performed on the FlOl swirl cup and splashplate assembly featuring the 
standard low pressure fuel injector, and the high pressure atomizing fuel 
nozzle. Both configurations were featured in some of the experimental test 
configurations. In each case care was taken so that the axisymmetric parts 
were concentric about the swirl cup axis. The splashplate assembly is asym- 
metric. The data obtained from these measurements are presented in Figures 35 
anid 36. Graphic 3-D presentations of Che data are shown in Figures 37 and 38. 
It is evident from these data that the fuel distributions are not axisymmetric. 
In both cases the greatest concentration of fuel lies in the region from 180 
td 360* clockwise aft looking forward. These results are in agreement with 
past results obtained on this swirl cup design. They point out the need for 
considerable detail in modeling the dome inlet if accuracy is to be achieved 
with analytical calculations. Some additional patternation data is presented 
in Figures 39 through 44. These are results from a similar swirl cup design 
ilJlustrating the effects of eccentricity among the axisymmetric parts, and 
of the fuel nozzle tip immersion through the face of the primary swirler. 
Neither of these effects were present in the experimental rig testing con- 
ducted as part of this program. However in a real engine, each swirl cup 
around the annulus can suffer from different degrees of eccentricity and fuel 
injector tip immersion. The resulting circumferential nonuniformity can have 
a significant impact on the measured combustor liner and turbine inlet gas 
temperature patterns. 

Airflow velocity measurements at the discharge plane of the FlOl swirl 
cup assembly were made using a five-hole (3-dimensional), wedge probe. An 
illustration of this probe is shown in Figure 45. In this test stand the 
swirl cup test piece is mounted onto a plate at the end of a 10.2-cm (4-inch) 
air supply pipe section. Air is supplied to the swirl cup through this pipe 
section at ambient temperature and atmospheric pressure. The wedge probe is 
positioned using a device which allows axial, vertical, and horizontal move- 
ment as well as rotation about the axis of the probe. Pressures from all five 
probe orifices are read on a bank of water manometers. The axial, radial, and 
tangential velocity components are calculated from these pressure measurements 
using a set of calibration curves generated for the five-hole wedge probe. 

The results of the wedge probe measurements made on the FlOl swirl cup assembly 
are shown in Figure 46. The probe was initially positioned at the swirl cup 
centerline at a distance of 0.25-cm (0.10-inch) downstream of the discharge 
plane. Measurements were made by moving the probe horizontally along the cup 
diameter away from the cup center at 0.13— cm (0.05— inch) increments. As^ 
observed from the results, at this measurement plane the swirl cup exhibits a 
recirculation zone extending radially outward 1.14— cm (0.45— inch) from the 
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Figure 40. Earlier Patternation Test Results Showing Effects of Eccentricity and Fuel 
Injector Tip Inunersion. 
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Figure 41, Earlier Patternation Test Results Showing Effects of Eccehtricity and Fuel 
Injector Tip Immersion. 
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Figure 45. 5 Hole (3-D) Wedge type Velocity Probe 
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center. The low tangential velocities reflect the offsetting effects of the 
counterrotation feature of this swirl cup design. At the center some small 
radial and tangential velocities were measured. The fact that these are not 
exactly zero may indicate that the probe was not exactly centered, and, or 
interference to the flow field was caused by the presence of the probe body. 
Although this type of intrusive measurement suffers from inaccuracies caused 
by the probe body interference, it is considerably faster and cheaper to 
perform than more sophisticated laser velocimeter techniques. The results of 
these measurements were easily modeled into the 3-D elliptic detailed internal 
flow model through the use of special swirler inputs incorporated in the 
coding. 
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6.0 FORMAL ANALYSIS OF INTERNAL FLOW MODEL CONCEPTS 


6.1 INTRODUCTION 

Comprehensive three-dimensional computer codes have been developed for 
the analysis of gas turbine combustors. The code prepared by the Garrett 
Turbine Company contains models for each of the relevant aspects of combustor 
flows. Northern Research and Engineering Corporation has assembled a model 
very similar to the Garrett model. These two models are incorporated in 
General Electric's internal flow module, INTFLOW. Another similar model 
named PACE developed at Rolls Royce has often been cited in the literature 
but is not available at General Electric. 

These programs represent major accomplishments in establishing the feasi- 
bility of computer programs which, at least in an approximate sense, include 
and codify physical submodels for a description of a practical gas turbine com- 
bustor. Many of the submodels, however, are inadequate and require major 
rethinking of both the relevant physics and the numerical algorithms. Hence, 
the simulation of the axisymmetric benchmark flows performed in this study uti- 
lized another code developed at General Electric, a versatile two-dimensional 
research code GETREF (general Electric T[urbulent Reacting Bjlow), which is dis- 
cussed in detail in this section of the report. 

Considerable efforts over the last 3 years have been made at the General 
Electric Company to improve the mathematical techniques and the physical sub- 
models. This inhouse work, demonstrated in this section, starts with a care- 
ful formulation of the overall conservation equations employing density- 
weighted averaging to obtain the mean equations. The appearance of additional 
terms in these equations compared to constant density flows is discussed. 
Alternative combustion submodels were developed accounting for either finite 
rate chemistry or equilibrium chemistry first with fluctuations of the chemi- 
cal source terms ("unmixedness") neglected. An equilibrium model accounting 
for these fluctuations by using an assumed probability density function (pdf) 
was developed and compared to the earlier submodels. Furthermore, improved 
turbulence submodels which account for anisotropies due to streamline curvature 
and swirl have been developed. The influence of density variations on the tur- 
bulence submodels was also considered. Means to reduce the numerical trunca- 
tion errors, which influence (to some degree) all finite difference solutions, 
have been incorporated into the computer codes. The proposed methods have been 
included in programs which treat only axisymmetric geometries, since the major- 
ity of the important physics can be incorporated and developed in 2-D computer 
analyses with a diminished complexity and cost compared to fully 3-D codes. 
Furthermore, most of the available benchmark quality flows are axisymmetric. 

The current generic mathematical techniques and physical models used for 
computation of turbulent combusting flows are analyzed and/or discussed in this 
section, while the application of GETREF to the selected axisymmetric benchmark 
flows is presented in Section 7.1. 
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6.2 CONSERVATION EQUATIONS 

The general equations which describe turbulent reacting flows are given 
below using Cartesian tensor notation for brevity. These equations are for 
the instantaneous flow field, Libby and Williams, 1980^^. (Repeated indices 
imply summation over that index. ) 


Overall Mass Conservation 




(13) 


Conservation of Momentum 


■it Ph-^-nr-’a 

J 1 J 
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Conservation of Energy (N Species) 
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Species Conservation 
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The mixture enthalpy h is related to the temperature in terms of the 
species enthalpies; 


N 

■■ ■ I ^ “i 

iii ‘ ^ 


(17a) 


‘i “ + h? ,hj - heat of formation of species i 


(17b) 


The equation of state for an ideal gas is given by 

p 

P 


»„T J (T^/ll.) 


(18) 


i=l 


Further discussion of the enthalpy and species conservation equations will 
follow. 

In combustor flows of interest, the term 3p/at in the enthalpy equation 

IS neglected on the basis that the Mach numbers in gas turbine combustors are 
low. 

In turbulent combustion systems, the time and length scales vary over such 
a wide range and the exact solutions (if available) would be so sensitive to 
initial and boundary conditions that it is not feasible to seek an exact ana- 
lytical or even a direct numerical solution of the above equations. For the 
high Reynolds number of interest in practical systems, the scales of the fluc- 
tuations are several orders of magnitude greater than those associated with 
molecular processes (diffusion). Turbulence reduces the scales of the fluc- 
tuations ultimately to a microscale range where molecular diffusion (viscous 
dissipation) becomes important, but this is an inertially controlled process. 
The models which follow, therefore, do not account for molecular effects, but 
rather assume the turbulence is inertially dominated. 

• standard approach to this problem is to decompose each independent 
variable (u^, h, Yj, etc.) into a mean plus fluctuating part, e.g.. 
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Tliis decomposition is carried out for each variable and then substitution is 
made into the governing equationsj with the resulting equations averaged in an 
ensemble sense (equivalent to time averaging for stationary flows). The aver- 
aging process, especially because of the nonlinear convective terms, and the 
chemical source terms in each species conservation equation, results in loss 
of detailed information and closure problems such as u'^ u'j are introduced, 
Kinze, 1959^^. Determination of such correlations is the essence of turbu- 
lence modeling. If conventional time averaging is adopted for variable den- 
sity flows , correlations (including the Reynolds stresses) such as p’u'j and 
p' u'i u'j arise. At the present time, insufficient data is available to 
accurately evaluate such terms. Instead, if density-weighted, or Favre-aver- 
aging, is utilized, viz. 


0 , '( 20 ) 


major simplifications in the form of the equations results. For example, the 
convection flux term in the momentum equations may be written 

'"‘I'j ' Vj 




n + a' 
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T— >* 




p n dt tad - 


P + P P'«'i * P'w'j + P'®'i«'j ( 21 ) 


» p iSjSj + a"£‘a"j] 

Obviously, the use of Favre-averaging simplifies the resulting expressions. 

The choice of either time-averaging or density-weighted averaging is not 
simply a matter of choice. Many investigators who have used time-averaging 
have neglected correlations including e' due to insufficient information. If 
that is done, or if Favre-averaging is used, equations arise which superfi- 
cially look like constant density equations where ( ~ ) correlations are 
replaced by ( ~ ) correlations. Additional terms areise in the modeled turbu- 
lence equations when correctly Favre-averaged that do not arise if e' correla- 
tions are neglected. 

I Favre-averaging provides a consistent framework to consider variable den- 
sity effects. Firstly, the density weighting results in equations describing 
the variation of the mean values of the conserved quantities Jones and Whi- 
tel aw, 1981^5.^ j[n the momentum equation, what arises is the mean momentum 
pui = = P u^. Secondly, as will be discussed below, the turbulence 

modeling of the Reynolds stress u"^ u"j at present can most directly be 
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carried over from constant density models which have been most severely tested. 
In essence, this is done with the "hope" that no significant errors result, 
Jones, 1979^®, despite the evidence of the necessity for density ratio 
effects in turbulence models, Libby and Bray, 1977^°, for premixed flames. 
Lastly, if density-weighted averages are compared with experiments, one must 
take care how the measurement_s_ were_taken before comparing data with calcula- 
tions, especially noting ^ - (p ^ lp’<p’/p, where <j) is any variable of Interest. 
Thus if conventibnal and Favre“averaged quantities are significantly different, 
the density fluctuation correlations are clearly important. 

The Favre-averaged form of the Equations (13), (14), (15), and (16) are 
given (for steady flows) by 


a 






( 22 ) 



^ (p S. 9 ) «-.j|- (p a-.?-") + (?) 


where ij) is any scalar variable (h, Y^, i - 1, . . . , N) and is the related 
soilrce term. Note that the molecular diffusion terras are ignored herein, 
under the assumption that the Reynolds number is sufficiently high. Also in 
any modeling it is assumed that Sc = Pr (in Equation (11), i.e., the Lewis num- 
ber is unity.) 

The means of handling the unknown velocity correlations pu "i u-j (i.e., 
the Favre-averaged Reynolds stresses), the velocity-scalar correlations pu"i (j)*', 
and the chemical source terms are discussed in detail in the following sections 
on turbulence and combustion modeling. The equation set (22) - (24) represents 
the mean flow equations of interest. 


6.3 NUMERICAI. METHODS 

This section discusses the basic features of the numerical algorithms 
used to solve the elliptic equations required to describe recirculating com- 
bustor flows; discretization of the convection term is emphasized. The need 
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solvers arises since real gas turbine combustors employ flame- 
Cabilizing recLrculacion zones including swirl flows as well or Irror rv 
reS^n"?;:""' ^ilufion air jers. In parfic\lar/ge„™:i%%™‘L foj ^fdacUr' 
sSwn numerical diffusion effects are discussed in detail, and results are 
shown to demonstrate the improvements gained. 

6.3.1 Current Method - Hybrid Scheme 

tic algorithm commonly used is based on the TEACH ellip- 

1C finite difference program, Gosraan and Pun, 1974^7 He.vplnr.of> at- t • 
Collage London. This method is well documentid i^the opei lUe?«Jrl®"a"d h.s 

198051. introduction to the method is given in Patankar, 

turbulence and turbulence/chemistry interaction modeling aspects will 
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inite volme, as opposed to strictly finite difference, approach. The orea- 

thr„OMnLrequftioLrt‘o®d«^rrjrLsure correctif 

si^LE (s^mi-impiicit Lthod for ppes^':;:-^Lh:d'” igokS^.^ai:;;;:? 

be soi:aV“^:r':t'iTe“ug\“%o^r““!y1':-;r'S:^ 

The procedure works iterativelv ba< 5 oH nn rt 

dent variablo Tho ranmom eased on initial guesses for each depen- 

aent variable. The conservation equation for each variable is successivelv 
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K H or Iteration. This form of linearization is stable, but 
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requires many line sweeps and iteration cycles to reach convergence. Conver- 
gence is usually judged to be reached when the sum over the entire computa- 
tional field of the absolute values of the residual error at each computation 
cell is reduced below some prescribed limit. For example, the sum total of 
mass continuity residuals is generally required to be less than 0.001 times 
the total inlet mass flux. The other residuals are also monitored. In 3-D 
flows some authors have targeted for a mass continuity residual of 0.01. 

Bpundary conditions are satisfied using thewall function method, Launder 
and Spalding, 1972^^. The viscous boundary conditions are satisified using the 
no slip condition at the wall and the law-of-the-wall to calculate the wall 
shear stress. 

The source terms are linearized in the finite difference form for the 
hydrodynamic variables according to S<j, = Sy - Sp^, where Sp > 0, Sp and S^^ can 
be functions of the coordinate and any dependent variable (based on a prior 
iteration). The source terms for tlie species conservation equations cannot be 
so linearized and, therefore, pose special problems and treatments. This topic 
will be treated more fully in the following section on chemical kinetics. 

The finite difference approximation employs so-called hybrid differenc- 
ing, i.e., second order central differencing is used for cell Peclet numbers 
(absolute value), e|u|Ax/u^, less than two; otherwise first order upwind dif- 
ferencing is applied. This methodology is required since second order central 
differencing of the convective terms for cell Peclet number greater than two 
results in oscillatory solutions, and often in divergence of the iterations if 
, a standard Gauss-Seidel procedure is applied. The result of this unwinding is 
that a numerical viscosity with magnitude e|u|Ax/2 is added to the turbulent 
eddy viscosity uj., resulting in excessive artificial or numerical diffusion 
in the flow field. This diffusional effect may exceed the physical turbulent 
diffusion and so plague the accuracy of the numerical simulations. For example, 
in the pioneering work by Allen and Southwell, 195553^ the first order upwind 
approximation was employed in a vorticity-stream function solution to the 
Kavier-Stokes equations for viscous flow around a cylinder. The results show 
that the downstream eddies are too short and they vary only slightly between 
R = 100 and 1000 on a coarse mesh. 

Furthermore, de Vahl Davis and Mallinson, 197654 pointed out that when 
the velocity vector is more than marginally skewed relative to the computational 
grid lines, and the cell Peclet number significantly exceeds unity in regions 
where diffusive transport normal to the flow direction is important, the error 
due to numerical diffusion may become so dominant as to obscure the effects of 
physical diffusivity on the flow. Also, Castro, 197955 showed that in regions 
of complex flow, particularly near sharp corners, the size of the truncation 
error associated with the first order upwind scheme cannot easily be reduced to 
insignificance, but it can be simply convected downstream leading to bad predic- 
tions over the rest of the flow field. In addition, Ralthby, 1976^® and 
Leonard 197957 have shown that in the presence of source terms, large errors 
may result from the first order upwind solutions for a convection dominated 
flow. 




Most of the published calculations for elliptic flows demonstrate various 
levels of agreement with experimental data. Often authors cite the need for a 
more sophisticated turbulence model to improve the quality of the calculation. 
However, in most instances, the amount of numerical error has not been ascer- 
tained, It is necessary to determine, for example, whether numerical diffusion 
has "smeared" the results in regions of high gradient. McGuirk, Taylor and 
Whitelaw, 198158 have published results for the calculation of the axisym- 
metric flow behind a bluff plate within a circular pipe. Grid refinement was 
carried out by examining for a given mesh the error due to false diffusion and 
then refining the mesh until such error was minimized. This is a laborious 
task, but one that should in some sense be automated to quantify the numerical 
errors involved. With three-dimensional calculations, the ability to demon- 
strate either grid independence of solutions or their degree of error, as sug- 
gested above, is severely limited. Coupled with the additional equations 
required for combusting flows, numerical testing for errors will be expensive 
and difficult, if not impossible, for three dimensional flows, Jones and 
Whitelaw, 1981^5. without combustion, some of these difficulties have been 
demonstrated by Green and Whitelaw, 1980^ and Green, 198l59. 

The difficulty of a numerical simulation is further compounded by the fact 
that for the cell Reynolds Number (Re = global Reynolds number, Ax = nondimen- 
sional mesh spacing) greater than unity, difference schemes of higher order 
formal accuracy do not necessarily promise smaller total error. Chen and 
Shubin, 1978^0 studied the one-dimensional, steady-state Burgers equation. 

They found that the error in computational results with formally second order 
accurate algorithms and coarse meshes varies widely, and does not increase as 
Ax 2 or (ReAx)2. Furthermore, the first order accurate algorithm can provide 
essentially the same solutions as do some of the second order algorithms; 
such results can be either better or worse than those offered by other formally 
second order schemes. Thus, they concluded that the formal order of accuracy 
of a difference algorithm may not reflect the magnitude of computational errors 
for large ReAx. Stubley, et al., 1980^1 demonstrated that it is not just the 
error introduced in the approximation of local function values or gradients by 
a particular discretization scheme which is important, but also the nature of 
the way this local error is propagated by the discretized version of the con- 
vection and diffusion terras in the differential equation. 


6 ,3.2 Alternative Methods for Three-Dimensional Flows 

Most of the numerical simulations done for the axisymmetric benchmark 
experiments employ the hybrid scheme since for two-dimensional flows the grid 
refinements, which make the cell Peclet always less than two (so that the sec- 
ond order central differencing is applicable), are possible within current com- 
puter capacities. For the three-dimensional flows of practical interest such 
grid densification in each coordinate direction is currently impractical so 
that the second order accuracy of central differencing cannot be realized. 
Stability considerations dictate the use of first order upwinding; however, 
this is now seen to be too diffusive to be of practical use. Clearly, a sec- 
ond order accurate scheme with an expanded stable range of Peclet numbers is 


106 





needed. This has been studied at General Electric, Shyy, 1983^2 gnd the 
results are presented below. Generally, all the alternatives to the hybrid 
scheme are necessarily more complex. Several potentially useful schemes are 
currently available. One such method, proposed by Rafthby, 1976^2^ is termed 
skew upwind differencing. The skew upwind differencing scheme, although like 
the conventional upwind scheme formally only first-border accurate, yields a 
significant reduction in numerical diffusion by taking the direction of the 
velocity vector into account. Another method, proposed by Leonard, 197957 is 
called QUICK (Quadratic Upstream Interpolation for Convective Kinematics). 

The QUICK scheme is based on the use Of upstream-shifted parabolic interpola- 
tion for every control volume surface on the computational grid and is free 
from (second order) numerical diffusion. Both skew upwind and QUICK have been 
examined by Leschziner, 1980^^, and by Leschziner and Rodi, 1981^5 for some 
idealized cases as well as for unconfined recirculating flows. The schemes are 
found to be superior to the conventional upwind simulation for the cases stud- 
ied. Both formulations have been shown, however, to involve boundedness prob- 
lems by under- and over-shoots in the solutions. It should also be noted that 
the skew upwind scheme does not resolve the difficulties associated with the 
source term when applying the first order upwind approximations; Leonard, 

197957 has shown that QUICK does give good resolution in the presence of the 
source term. On the other hand, however, if an iterative method such as Gauss- 
Seidel is adopted for the solution procedure, Han, et al. , 1981®° reported 
that special care must be taken for the QUICK scheme to guarantee convergence. 

Another finite difference approximation which has attracted comparatively 
little attention was the one used by Atias, et al., 1977^7, They studied a 
second order upwind scheme for discretizing the convection terms in the vortic- 
ity transport equation and found that it has the potential of yielding suffi- 
cient accuracy. No systematic study has been made for this scheme as yet; but 
Gupta and Manohar, 1978^® poihted out the restriction of numerical stability 
imposed by the von Neumann type of analysis. 

In the study conducted at General Electric, Shyy, 1983®2j attempts have 
been made to clarify some of these ambiguities associated with the numerical 
simulation of a high cell Peclet number flow as well as to find a comparatively 
accurate finite difference scheme and an effective procedure to get the solu- 
tion to the discretized equation. The restraints set by the size of the cell 
Peclet number on the effectiveness of a given finite difference scheme are 
studied at first. The accuracy of a finite difference approximation to an 
idealized flow problem in the presence of a source term is then analyzed. Com- 
parisons are made by studying the solutions by different schemes for several 
test problems. The effects of numerical diffusion due to the inclination of 
the computational grid lines in a multi-dimensional flow field are also com- 
pared among different schemes. The influences of different boundary condi- 
tions on the mathematical accuracy and on the physical reality of the numeri- 
cal solutions are examined. One objective here is to shed some light on the 
relative merits of several schemes and to contrast their performances under 
both "mild" and "stringent" flow conditions. Five different finite difference 
schemes approximating the convection terras are analyzed: the first order 

upwind, the skew upwind, the second order central differencing, the second 
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order upwind, and the QUICK' scheme. Test cases are chosen to be simple but 
relevant. A stable iterative procedure is constructed to give convergent solu- 
tions for all the schemes considered. 

In what follows, exclusive use of the second order central differencing 
will be adopted to approximate the diffusion terras in the governing equation, 
e.g., 




i.j 


A,2 


+ T, 


(26) 


where is the truncation error resulting from replacing the diffusion term 
with the finite difference approximation; the notation used in Equation (26) is 
shown in Figure 47. 

For the finite difference approximation to the convection terras, the fol- 
lowing schemes will be considered and tested: 

a) first order upwind 


”i>j *i,i ~ 


I 

dx i.j 


Az 


+ T^.for n > 0 




(27) 


where is the . truncation error inherent in replacing the convection terra 
with the finite difference approximation. 

b) second order central differencing 


a(it») I _ °i I.j + l,j 
?z i.j 2Az 


(28) 
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c) second order upwind 
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3M) 


I . 


nr ♦ Tc f»t • > 0 


(29) 


i> j 


d) QUICK 


lll<-V2.j'i*2,j * * T. '"t ■ < « 


d(nf) 


* K-2.j»l-2.j) * t” » > 0 


15“ 'i.j 


-1- ( -J:« e + ^ 4 - ^ f 




(30) 


Skew upwind differencing is similar to the original first order upwind scheme 
for a one-dimensional problem. For a multidimensional flow problem, it can 
reduce numerical diffusion by taking the direction of the velocity vector 
into account as will be described later in this report. 


6.3.3 Numerical Accuracy and Cell Peclet Number 

It is well known by now that numerical schemes (to any finite degree of 
accuracy) introduce numerical diffusion and dispersion in roughly the same way 
as physical diffusion and dispersion in the phenomena of fluid flow, Richtmyer 
and Morton, 1967^^, Chu, 1978'®. This can be examined by expanding the finite 
difference equation in a Taylor series, to get the original differential equa- 
tion plus higher order terms which represent the truncation errors introduced 
in the course of approximation. The resulting equation is called the modified 
equation. Table IX shows the coefficients of the first six derivatives in the 
modified equation 
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q# + 

' XXX 


d# 

XX 


ed 


(4) 


+ 


( 6 ) 


+ HOT 


(31) 


where HOT represents the higher order terms in Taylor series, derived by the 
combination of one of the finite difference approximations discussed in Eqs. 
127) to (30) to the convection term and Eq. (26) to the diffusion term of the 
one-dimensional , steady state linear Burgers equation (source free) 


“ ^^xx' ®» \) * coaxtant* > 0 


(32) 


If we expand the dependent variable $ in Eq. 
the x-direction we obtain 


(31) as a Fourier series in 


r ik X 

2 e J 

J*— • J 


(33) 


_ The Fourier component of the shortest wavelength resolved by a finite 
difference mesh is of wavelength 1 = 2ax. The corresponding wave number is 

wavelength is 1 ^ 3 ^= L, which is the total length 
panned by the mesh. The corresponding minimum wave number is = 2 tt/L. 

If 'one substitutes a Fourier component (pj = a^e^ into Eq. (13), then one 
can easily find that all the schemes in Table-'lX except the first order 

" ^in with the errors pro- 

viel^reL^LM " ' On the Other hand, as k = no scheme can 

yield reasonably accurate approximations. This can be seen from Table X which 

shows the values of o and 3 for the Fourier components with k “ k,„av for the 
various schemes, where "max* tor cne 


o = 


q^ + 

^ XXX *^^x 


(34a) 


(d-\))dj^ + 


\)P 


XX 


(34b) 
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Table X. Values of a and 3 in Eq. (34) by Substituting 
the Fourier Component of Wave Number k * tt/Ax 




All the values of a and 0 are at least of the order of imityj in the 
region where those high wave number Fourier components are important for the 
solution, the numerical diffusion and dispersion are by no means negligible 
even though all the schemes studied here, except the first order upwind, make 
exact approximations to the first and the second derivatives in Eq. (32). 
Furthermore, this difficulty cannot be resolved by using a formally higher 
order scheme. 

For a high Peclet number flow problem, the convection terms are dominant 
in the main part of the flow domain; the real solution might vary rapidly in 
some thin regions where the convection terms are balanced by the diffusion 
terms due to, for example, the restraint imposed by the boundary conditions. 
Hence, in a numerical viscous flow simulation the mesh size should be fine 
enough so that the ratio of ucfjjj/ (V(j)j^) in Eq. (32) is of the order of unity, 
at least at the highest resolvable wave number, k = k-^^ax “ r/Ax. If the ratio 
is much larger than one, then the convection term can never be balanced by the 
viscous term in all the numerically resolvable scales. If a Fourier component 
(j)j = a^e^^j^ is substituted into u({>x/ (^<(>xx) “ Knav used, then as 

pointed out by MacCormack and Lomax, 197971, the condition for a unity ratio 
between the convection and the diffusion terms is 


eAz 


Jft- 0(1) 


(35) 


At high Peclet numbers, the mesh size cannot, with present computers, be 
made fine enough to fulfill this condition. The viscous diffusion at wave 
numbers higher than those resolved by the numerical approximation must be 
accounted for through modeling. Hence, a suitable finite difference scheme 
for a high Peclet number flow should be the one that maintains good accuracy 
in the convection dominated region; in the regions where the convection terms 
and the viscous terms in the original differential equation must balance each 
other, but numerically this cannot be done exactly, the truncation error should 
enhance the weight of the viscous terras. By doing this, though the fine struc- 
ture of the real solutions cannot be resolved due to the large error introduced 
in high wave number component part, the results obtained from the numerical 
apprcDximation are, in general, physically correct and useful; the disturbances 
arising from an inaccurate numerical simulation in a thin layer can be damped 
out without propagating into the main region. As will be shown later, the 
truncation error of a given finite difference scheme may or may not be able to 
increase the viscous effects where needed. Hence, a formally higher order 
scheme does not necessarily perform better for the high cell Peclet number 
problems, due to the inadequate resolution of the rapidly varying solution 
in the thin layers. 


6.3.4 Test Problem I: Bou nd ary Layer Type Flow 

Eq. (32) with the boundary conditions = 0 at x = 0, and (|> = 1 at x = 1 
is considered here. The exact solution of this problem is 
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#(x) 


1-e 


Pe* 


1-e 


(36) 


where - ul/v is a Peclet number. For small to moderate pg, (J>(x) displays 

smoothly throughout the entire domain. However, as is 
increased much larger beyond unity, the solution becomes one of a boundary 
layer type in that 4>(x) is virtually zero except in the region near x = 1 
within a thin layer of thickness ?= 1/p^ wherein the entire variation in the 
solution is contained. Upstream of this boundary layer, the flow is entirely 
convection dominated. The numerical solutions to this test problem by the 

schemes described in Eq. (27) to Eq . 

UO) for different values of cell Peclet number are contained in Tables XI 
through XIII. For PeAx.= uAx/v =0.2 the soulutions by QUICK and the second 
order central differencing are the most accurate while that by the second order 
upwind IS also very acceptable. For PeAx less than 1, the Fourier components 

T. T r *^ax are not crucial for ^ in 

Eq. U6), and the leading truncation error term of a finite difference approx- 
imation IS representative of the order of total numerical error. In Table XI 
the first order upwind yields the solution with the largest error due to the 
introduction of the numerical viscosity (see Table IX) to the second derivative 
term in (^2). It is noted that among the four different schemes described 
in Eqs. (27) to (30), only the first order upwind always gives a diagonally dom- 
matrix for the set of difference equations. For p^Av = 10 
and 100, both QUICK and the second order central differencing show solutions 
mth wiggles; the magnitudes of this spurious oscillation in the solution by 
central differencing are generally much larger than those by 
I ’ throughout the whole range of p^^^, the solutions by both the 
first order upwind and the second order upwind are wiggle free. This can be 
analyzed by studying the roots of the characteristic equation associated with 
a specific scheme. For example, if QUICK is used, the resulting finite differ- 
ence equation for Eq. (32) at the grid point i is ^ 


eAx 


i - 2 


“ 8 ^i-1 - (1 - P 


eAx^**i+l ~ ® 


(37) 


The exact solution o|f Eq. (37) is given by Dahlquist and Bjorck, 197472 


'i ■ * T2^r' * Vs 


i-1 


(38) 
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Table XI. Numerical Solutions J;p Test Problem I. 
(Number of Grid Points = N + 1 = 11, = 0.2) 


Metbod 


*1 ^ ^ ^ U U ^ U ^ ^0 


Eq. (26) 


0 3.47 7.70 0.129 0.192 0.269 0.363 0.478 0.619 0.790 1.00 

E-2 E-2 


Pint order 0 3.85 8.4 8 0.140 0.207 0.287 0.383 0.4 98 0.6 3 6 0.801 1.00 

npwind E-2 E-2 


Second order 0 3.58 7.88 0.131 0.195 0.272 0.367 0.482 0.622 0.792 1.00 qq 

npwind E-2 E-2 S SB 

: 


Second order 0 3.45 7.67 0.128 0.191 0.268 0.362 0.477 0.618 0.790 1.00 

central E-2 E-2 

OPICE 0 3.53 7.77 0.129 0.193 0.270 0.364 0.479 0.619 0.791 1.00 

»» — 


<0 3 
c > 
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Table XII. Numerical Solutions to Test Problem I (11 Grid Points, = 10.0) 


Method 

*1 

^2 

*3 


*5 

*6 

*1 

^8 

^9 

^10 

Eq. (26) 

0 

8.19 

E-40 

1.80 

E-35 

3.98 

E-31 

8.76 

E-27 

1.93 
E-2 2 

4.25 
E-1 8 

9.36 
E-1 4 

2.06 

E-9 

4.54 

E-5 

First order 
upwind 

0 

3.90 

E-10 

4.65 

E-9 

5.12 

E-8 

5.64 

E-7 

6.21 

E-6 

6.83 

E-5 

7.51 

E-4 

8.26 

E-3 

9.09 

E-2 

Second order 
upwind 

0 

-3.46 

E-10 

-3.57 

E-10 

3.34 

E-9 

6.59 

E-8 

1.05 

E-6 

1.65 

E-5 

2.59 

E-4 

4.07 

E-3 

6.38 

E-2 

Second order 
central 

0 

-4.41 

E-2 

2.21 

E-2 

-7.72 

E-2 

7.17 

E-2 

-1.52 

E-1 

1.83 

E-1 

-3.19 

E-1 

4.35 

E-1 

-6.96 

E-1 

QUICK 

0 

1.06 

E-5 

1.27 

E-4 

-2.28 

E-4 

9.22 

E-4 

-2.79 

E-3 

9.21 

E-3 

-2.96 

E-2 

9.58 

E-2 

-3.09 

E-1 
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Table XIII.. NuroericaL Solutions- .to_ Test Problem I (11 Grid Points, = 100) 


Method 


Eq. (26) 


First order 
opvind 


Second order 
n]nrind 


Second order 
centrsl 


QUICK 




I 

0 

B 

I 


-5.05 

E-11 


-4.15 


u 

% 

u 


*1 



0 

0 

0 

0 

0 

0 

0 

1.01 

E-15 

3.19 

E-15 

1.45 

E-12 

9.69 

E-11 

9.61 

E-9 

9.71 

E-7 

9.80 

E-5 

-9.06 

E-11 

-1.40 

E-10 

-2.11 

E-10 

-3.03 

E-10 

1.47 

E-9 

2.92 

E-7 

4.41 

E-5 

1.69 

E-1 

-4.32 

3.53 

E-1 

-4.52 

5.52 

E-1 

-4.72 

7.67 

E-1 

2.51 

E-3 

-2.65 

E-3 

8.94 

E-3 

-1.70 

E-2 

4.09 

E-2 

-8.85 

E-2 

2.01 

E-1 


0 1.00 
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where Zi, Z2, and Z3 are zeros of the characteristic equation 

^ * kA.' ^ * 4 P.AX>Z ^-<1 - |p.A.>^ 


These zeros are 


( 39 ) 


1 , 







(1 - 


T P.Ax> 


he^constants yi, Y2, and Y3 are determined by the two boundary conditions 

- 1. and the value of ^2 which is calculated by a starting calcu- 
lation method. Here, the first order upwind may be used for both QUICK and the 
second order upwind schemes; it is found that, for this specific problem, th^ 
solutions obtained in this way are comparable with other more accurate strate- 
gies such as the one that will be discussed later on for the flow problem with 

indicates that, for the problem considered, QUICK 

in Eq. ( 40 ) has a negative sign for p^av > 8 / 3 . Hence d a = fi/'i io i-u 

critical cell Peclet number for QUICK. Table XIV shows the roots of the char-^ 

ciatld^wi^h^Jr*^^^"- critical cell Peclet number asso- 

thP n J various schemes for the test problem I. Those roots of both 

the first order upwind and the second order upwind are always positive- the 
solutions for the model problem by the two schemes are wiggle free. The criti- 
cal cell Peclet number for second order central differencing is, as is well 
As PeAx IS far longer than unity the two characteristic roots of 
e second order central differencing are very close to each other in magni- 
tude. Hence, as was shown by Gresho and Lee, 197973 , the magnitudes of the 

= dependent upon whether the total number of mesh 

f Nevertheless, it is clear that, for this model problem, 
wi-n PeAx. both the second order central differencing and QUICK 

11 not be very satisfactory; the higher order terms in the modified equations 

J!:: adequately to balance ?he co°n- 

™ ^ boundary layer without letting the disturbances propagate 

nto the main region, therefore, the numerical solutions oscillate. On the 
other hand, the solutions by both the first order upwind and the second order 
upwind are wiple free and acceptable: the second order upwind derives more 
accurate simulations than the first order one for the whole range of p * . 

bbe truncation errors in these two schemes are able to help damp out 
the disturbances where needed more effectively than those in QUICK and the sec- 
ond order central differencing. By examanlng the performance of the three sec- 
ond order accurate schemes, it is clear that the formal order of accuracy loses 
ips meaning if a boundary layer exists in a high PeAx problem since the leading 
truncation error term no longer represents the real size of the error intro- 
duced by the numerical approximation. But for a low PgAx flow, the formal 
order of accuracy is still a good basis by which to judge the relative accuracy 
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Scheae 


Boots 


Critical Value of 


•'eAi 


First Order 
Opwisd 


1 + p 


eAx 



Second Order 
Central 


1 + 


"eAx 


1 - 


"eAx 


QUICK 


( 1 ) ± 


( 1 


2 




between two, schemes with different order, e.gt , the first order upwind and 
QUICK, 

Another aspect concerning the validity of the numerical solutions with 
the large values of the cell Peclet number which is worth commenting on is 
that no scheme is accurate if a quantity such as ({ix is to be calculated at 
the right hand boundary point for this model problem. This is, apin, due to 
the fact that the high wave number Fourier components, which are important at 
right hand boundary for high Pe^x cannot be evaluated accurately in 

the course of the numerical approximations by a finite difference scheme. 

Here, (jiN+i » 1 and the smallest value of that one would like to see is 
zero (otherwise the wiggles start to appear). Hence, the numerical calcula- 
tion for (j)x at X = 1 is, at best, (<}>n+1 ” “ 1/Ax. This result is inde- 

pendent of the Peclet number, in contrast to that which was obtained from the 
analytical solution which increases directly with PeAx* 


6.3.5 Iterative Method 

The numerical procedures suitable to get the solution of the finite dif- 
ference equations approximating the lai’ge Pe^x flows are discussed. For test 
problem I, if the standard SOR procedure is applied, it is learned that a con- 
vergent solution for any large values of Pe^x 1*® guaranteed only for the 
first order upwind scheme; all the other three schemes fail to provide conver- 
gent solutions if Pe^x 1® larger than some values. The reason is as follows. 
For a system of real linear equations X<|) = in which the matrix A can be 
written as A=I-L-U, where L is the strictly lower triangular matrix, U is the 
strictly upper triangular matrix, and I is the identity matrix, the Jacobi 
matrix Gj and the SOR matrix G^, corresponding to the relaxation factor u, 
can be written as 


Gj - L + U 


(41a) 


G ■ (I - dL)~^ 1(1 - »)I + 
10 


(41b) 


If the matrix A is consistently ordered, then the following relation 
exists between the eigenvalues p of Gj and the eigenvalues \ of G^^, Varga, 
196274 

(X + » - 1)* ■ #^|i^X (^2) 


121 




ORIGINAL PAGE rs 
OF POOR QUALITY 


Young, 1971^^ showed that an iterative method is convergent if and onlv 
< “maxi where ^ 


■** 11 + Mj (1 - 


(43a) 


n - Hg + r^» 


“I 


(43b) 


Furthermore, if either -0 or pj » 0 and |pr!<|, then the following for- 
mula IS available to determine the optimum value of the relaxation factor 
“opt Varga, 1962'4 » 




(44) 


We now compare the eigenvalues of Gj from the second order central dif- 
ferencing scheme with those from the first order upwind scheme for the test 
problem I. The central differencing approximation for Eq. (32) at the noint 

1 IS ^ 


^ - T »i-i * <1 - 


(45) 


The first order upwind approximation for Eq. (32) at the point i is 




i ■ ■* W ^ - 1 ♦ *i . ll 


(46) 


The eigenvalules p of the Jacobi matrix for both schemes 


are: 


Central Differencing 


(m) 


central 


I 1 - (^^)^ cos , k » 1, . . . , N 


(47) 
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(M) 


1+p 


Fn 




kn 

2 N+1 


1 . 


N 


(48) 


Since all (u)jTu in Eq. (48) are always real and less than one in magni- 
tude, It is known from Eq. (43a) that any value of the relP'^ication factor rang 
ing from 0 to 2 is allowed to give the convergent solution. For the second 
order central differencing scheme, all (u)central (47) are real and 

less than one in magnitude for PgAx 2; they become purely imaginary for 
PeAx >2. Hence, one learns from Eq. (43a) that, for the second order cen- 
tral differencing scheme, is 2 for pg^x < 2 and is less than 2 

^o>^.PeAx 2. Furthermore, if pg^^ » 1. then Eq. (44) indicates that the 
optimum value of the relaxation factor is 


(u y • s (u , s i 

opt'c.ntr.1 ••*)„„„, (49) 

(H— J— ) 

which is a very small number. The eigenvalues of Gj from the second order 
upwind and from QUICK are not known analytically. But it is proposed here 
that a relaxation factor u, analogous to that in Eq. (49) of the form 



a 


(50) 


be used for the iterative process to get the solutions of the various schemes 
considered here, a is the sum of the absolute values of all off-diagonal ele- 
ments in the ith equation (normalized with respect to the absolute sum of the 
diagonal elements) if the i^h point is to be calculated. For example, the 
value of (1 +0:5 PeAx). (2 + 2.5 PeAx>/(2 + 1.5 PeAx) » and (2 + 11 PgAx/8)/ 
j PeAx/°7 IS used for the second order central differencing, the second 
order upwind, and QUICK, respectively in the test Problem I. The iterative 
method with those under-relaxation factors has produced convergent solutions 
by all the schemes for all the test problems studied here. Of course, for the 
first order upwind scheme, the standard SOR procedure is always effective in 
giving convergent solutions j SOR should still be used for the first order 
upwind. 


No attempt is made here to compare in detail the computing time needed to 
get the final solutions for those finite difference equations using the relax- 
ation factor just described. But for all the problems with high PeAx tested 
in his study, including the one-dimensional convection diffusion flow, the 
multidimensional flow, and the flow with a source term, the second order upwind 


123 


t ! 

‘ i 


scheme and QUICK require almost the same number of iterations which are 
slightly more than those needed for the first order upwind but are less than 
those for the second order central differencing. This is contrary to rhe 
findings of Atias, et al., 1977^^ who applied the Gauss-Seidel method to 
solve the two-dimensional vorticity transport equation and reported that the 
second order upwind is always somewhat less efficient than central differ- 
encing. For the problem they solved, whether the present iterative method 
will show the same relative efficiency between the two schemes or not is open 
to further investigation. Nevertheless, in the study by Atias, et al., upper 
limits exist for the cell Peclet number for both schemes using the Gauss-Seidel 
method, as just explained. Another aspect worth commenting on is that the 
numerical stability restraint (in von Neumann sense) of an explicit scheme for 
the time-dependent problem is not equivalent to that of an iterative scheme for 
the steady-state problem. Hence, the criticism by Gupta and Manohar, 197868 
regarding the stringent stability limit of the second order upwind scheme for 
the time-dependent equation is not applicable here. An illuminating compara- 
tive study for the various iterative methods applied to large problems with the 
convection terms approximated by the second order central differencing scheme 
was given by Botta and Veldman, 198276. 


6.3.6 Test Problem II; Flow with Source Term 
Leonard, 1979^7 chose the model problem 


N) ® constant* > 0 
^(0) « 0 

x-L - 0 


(51a) 

(51b) 

(51c) 


i 

y 


B 



to study the accuracy of the finite difference approximation to high PoAv flow H 

in the presence of a source term. Further investigation is proposed for the ^ ‘ 

same problem by varying the distribution of S(x) to set up the limit of an 
adequate numerical simulation as well as comparing the performances among dif- 
ferent schemes. The basic form of S(x) used here is f 
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S(x) » »x + b, 0 1 X j( Xj 


(52) 


as shown in Figure 48. 

For the extreme case of PeAx"*'"> Leonard, 1979^7 has shown that the solu- 
tions given by QUICK are much better than those by the first order upwind which 

follow closely the exact solution of Pe^x “ those by the second order 

central differencing which are more susceptible to the changes of the down- 
stream boundary condition; here we concentrate on the comparisons between the 
solutions by the second order upwind and by QUICK for PgAx ~ 10® • For a source 
spanning several mesh lengths, both schemes can give accurate solutions for 
Pe4x"*’“> only complicating issue is the strategy for the starting calcula- 
tion. Figure 49 shows the solutions calculated by both schemes using the first 

order upwind scheme to calculate the first unknown value of the dependent vari- 

able, ip2- Figure 50 shows the solutions by both schemes putting the left hand 
boundary point, x = 0, at the mid-point of the first numerical mesh, and assign- 
ing the value of <{i at the first grid point, <t>i, by linearly extrapolating the 
values of if 2 (which is located at x = 4x/2 and the boundary value at x = 0. 

Here, cj) = 0 at x = 0, hence cf) 2 ^ = -cp2‘ The solutions calculated in this way, as 
suggested by Leonard, 1979^7^ are more accurate than those shown in Figure 
49 for both schemes. Also, in both Figures 49 and 50, the accuracy of the solu- 
tions by the two schemes are comparable. Figure 51 shows what changes can be 
made for the numerical results if f = 0, not ~ 0, is used as the downstream 
boundary condition. For this example, the solution by the second order upwind 
scheme is more satisfactory than that by QUICK. The restraint imposed by the 
downstream boundary condition forces the QUICK solution to show oscillations 
because the higher order terms in the modified equation fail to help damp out 
the disturbances in that thin layer. Figure 52 compares the two numerical 
solutions with the exact solution for a more rapid varying source term spanning 
two mesh lengths. Neither of the two numerical solutions is totally accurate 
for all grid points but the one by second order upwind is better. Figure 53 
shows the numerical solutions by the two schemes against the exact solution 
for a more stringent source distribution. For this problem the Fourier compo- 
nents with the wave number equal to or higher than ir/4x are- important to the 
exact solution, and, as was discussed previously, tkey cannot be resolved sat- 
isfactorily by the finite difference approximations; the solutions by both 
schemes carry unacceptably large errors even though the exact value of <f2 
assigned to eliminate any error involved in the starting calculation method. 
This demonstrates the fact that any attempt to accurately simulate a flow with 
the length scale 2Ax (or less) is out of the question; a more refined mesh must 


(tx« + b) (x, Xj) 

= X + - — ; + b) , *- < X i Xj 

*2 *2 ^ ^ 


+ X., 
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Numerical Solutions to Test Problem II (Downstream 
Boundary Condition: d *5/d x - 0) Starting Calculation 
Method: (-AX/2) = -i>i ^2’ first order upwind. 
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Figure 51. Numerical Solutions to Test Problem II (Downsteam 
Boundary Condition: ^ = 0) Starting Calculation 

Method: = ^ (- Ax/2) = -^2> by First Upwind 
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Figure 52. Numerical Solutions to Test Problem II (Downsteam 

Boundary Condition: d ^/d x*0) Starting Calculation 
Method: First Order Upwind. 
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be used for all schemes considered here. On the other hand, for a problem with 
a tnild source distribution and downstream boundary condition such as that shown 
in Figure 50, an accurate numerical approximation is possible for large pgAx 

because the diffusion term in Eq. (51a) plays no important role in the whole 
flow region. 


^•3-7 Test Problem III; Numerical Diffusion Du e to Streamline-to-Grid 
Skewness ~— 

To study the problem of numerical diffusion in a multidimensional flow, 
the following idealized case is considered at first: 


* 0 , 0 i X 1 1, 0 X y 1 1 

(53a) 

9(0, y) - 100y“, 0 < y i 1 

(53b) 

f(x, 0) - i00(-5^)“, 0 1 X i 1 

(53c) 


where cot e = y/v, and u and v are positive constants. The difference between 
the first order upwind and the skew upwind schemes is analyzed here to show 
some basic characteristics of numerical diffusion, The first order upwind 
approximation to Eq. (53a) at the grid point (i, j) can be written as follows: 


where 


(54b) 


Eq. (54a) shows how is calculated by linearly interpolating between 4>i_i 
and <j)£ at the point where the straight line connecting the grid points 
(i-l, j) and (i, j-1) intersects the velocity vector at the point (i,j). It 
13 known that when convection dominates in establishing the spatial distribu- 
tion of the quantity <t>, the gradient in in the streamwise direction vanishes, 
and the cross-flow gradient depends on the upstream boundary conditions. From 
Eq. 154b) It IS clear that the first order upwind gives the exact solution if 
the exponent n in Eqs. (53b, c) is unity; otherwise, the error resulting from 
the linear interpolation in Eq. (54a) will cause numerical diffusion in the 
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cross stream direction. Recognizing this, Raithby, 197663 devised the fol- 
lowing scheme to approximate Eq. (53a) more closely: 


Ax Av 0 


where 










Eq. (55a) can be rewritten as follows:' 


'i.j ■ Vi.j-1 "iiS'Vi.j - 


Figure 54 compares the way to calculate c()ij by the first order upwind and 
that by the skew upwind. It can be seen that, although both schemes are for- 
mally exact only to a linear variation of ♦ profile in the cross-stream direc- 
tipn, the calculated by the skew upwind can be improved substantially 

upwind. As for the second order upwind 
and QUICK, both are exact to the quadratic profile of <D in the cross-stream 
direction. Figure 55 shows the numerical solutions to Eqs. (55) by the four 
schemes: first order upwind, skew upwind, second order upwind, and QUICK, foi 

m=2 and Cot 0-2. The skew upwind is used to calculate the first unknowns fror 
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4 in Upstream 
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of 
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Figure 54. Comparison of Accuracy Between First Order Upwind and Skew Upwind Approximations 
to Eq. (53): Lines 1 and 2 are Farrell* 
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the upstream boundaries for both the second order upwind and QUICK. In addi- 
tion, the downstream boundary conditions are also needed by QUICK, and the skew 
upwind is also used there. In Figure 55, the grid points connected by any sin- 
gle streamline should possess the same value of 4 ,. While the first order upwind 
gives the solutions with some noticeable numerical cross-stream diffusion, all 
other three schemes give very good approximations. It is noted that the small 
errors appearing in the solution by QUICK are caused by the calculation for the 
downstre^ boundary values. Figure 55 demonstrates that the first order skew 
upwind differencing can reduce the cross-stream numerical diffusion substan- 
Figure 56 shows the numerical solutions by the four schemes for the 
same equation but with a more stringent upstream boundary condition; a step 
function, instead of a n-th power polynomial, is used. Here, again, because of 
the importance of the high wave number Fourier components to the real solution, 
no scheme is able to give a very accurate approximation; the sharp variation in 
the real solution is smeared out due to the numerical diffusion. Also, the 
solutions by the skew upwind, by the second order upwind, and by QUICK all 
exhibit overshoots; the solution by QUICK also exhibits undershoots. Globally 
speaking, QUICK and the second order upwind perform better here. 

After examining the possible error caused by the numerical approximations 
to the convection terms, a two-dimensional flow problem including both the con- 
vection and the diffusion terras, as stated in the following, is considered next. 


- \)(^„ + #yy) .0 i X i 1. 01 y 1 1 (57a) 

9 ( 0 , y) * 100, 0 < y 1 1 ( 57 b) 

^(x, 0) ■ 0, 0 1x11 (57c) 


where u and v are positive constants; u/v = 2, and the cell Peclet number based 
on uAx/v = 40 are used in the study. Two different downstream conditions are 
investigated. The periodic values of 4 ) along a streamline, or the zero value 
of 4 >. Five numerical schemes, first order upwind, skew upwind, second order 
central differencing, second order upwind, and QUICK, are compared in this test 
problem. For the first of the downstream boundary condition, i.e., the periodic 
values along a streamline, a typical profile of (j, plotted as a function of y at 
X|= 1 -2ax, given by the numerical approximations are shown in Figure 57. The 
profile by the first order upwind is the smoothest one due to the excessive 
numerical diffusion. The solution by the second order central differencing, on 
the other hand, shows the wiggles in the flow region. The other three schemes 
all produce the solutionis with some limited extent of overshoots near the down- 
stream boundary. Among them, the solutions by QUICK and by the second order 


Figure 56. 


0 0 

Numerical Solutions to Eq. (53) (Upstream Boundary Conditions: 
0 (0,y) = 100 0< y< 1, 0 (x,0) = 0, 0< x< 1) . 










upwind are the rapre satisfactory. Compared with the true solution for the first 
kind of boundary condition, the true solution for the second kind of downstream 
boundary conditions, i.e., (ji = 0, and 'high flow should be of the charac- 

ter that a boundary layer region is formed in the vicinity of the downstream 
boundary v/ith the solution in the main region of the flow field virtually unaf- 
fected. The solutions by the five numerical schemes for this problem are com- 
pared in Figure 5b. Those by the second order central differencing and by 
QUICK show unacceptably large magnitudes of oscillation. The solutions by the 
three upwind schemes are, compared to those in Figure 57, unaffected by the 
change of the boundary condition except for the points far downstream. In this 
case, the best scheme is the second order upwind which shows a reasonable com- 
promise between an accurate simulation in the convection dominated region and 
an effective enhancement to the diffusion term in the boundary layer region 
which prevents the disturbances from propagating. 


6.3.8 Summary and Conclusions 

Five different finite difference schemes: first order upwind, skew 

upwind, second order upwind, second order central differencing, and QUICK, 
approximating the convection terms in the equation governing fluid motion have 
been studied. It is shown that in a finite difference simulation with large 
cell Peclet number, the high wave number Fourier components of the real solu- 
tions cannot be evaluated accurately. As a result of this, in a convection 
dominated flow, if the viscous terms are required to balance the convection 
terms in a thin layer close to the downstream boundary due to, for example, a 
Dirichlet type boundary condition being applied there, a finite difference 
approximation can, at best, use the truncation errors of the approximations to 
help damp out the disturbances in that thin layer. The detailed structure of 
the real solution is unresolved. Even this may or may not be accomplished by 
a given scheme, and hence for a high cell Peclet number flow the formal order 
of accuracy is a poor criterion to judge the performances among different 
schemes. Among the five schemes examined in this study, the second order 
upwind scheme gives the most . satisfactory results for the particular test 
cases in this study. It still, however, exhibits overshoots in the solution 
to a limited extent. On the other hand, although both the second order central 
differencing and QUICK are formally of the same order of accuracy as the second 
order upwind, they fail to enhance the viscous terms properly, in the region 
where needed, for a high cell Peclet number flow problem, and noticeable spu- 
rious oscillations appear in the numerical solutions. It should be noted that 
the magnitudes of those spurious oscillations in the solutions by QUICK are 
generally less serious than those by the second order central differencing. 
Furthermore, if no boundary layer type of region exists in the real solution, 
the accuracy of the approximating solutions given by QUICK and by second order 
upwind are comparable. As for the first order upwind, it is free from produc- 
ing the unphysical over- and under-shoots in the solutions for all the test 
problems, but fails to give accurate approximations in the presence of a source 
term and shows too much numerical diffusion in the convection dominated region 
for a multi-dimensional flow. The skew upwind is able to reduce numerical dif- 
fusion by the first order in the cross-stream direction substantially, but it 
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Figure 58. Numerical Solutions to Eq. (56) at x = 1 - 2A x 
(Downstream Boundary Conditions; 0-0 ), 
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skew upwind scheme also fails to to liy ^ second order upwind. The 

the source terms, Raithby, 197663 For*'airth^*^f the presence of 

rate schemes studied here, the standard SOR ,^?™ally second order accu- 
large values of p^a • a rnnHifioa i ^ technique cannot be applied for 
tions for all cases^considered “at ion factor yields convergent solu- 


6.4 TURBULENCE CLOSURE MOnRT.g 

^ This section discusses the eeneral 

points in the field of turbulencf modeling specific 

one point correlations, u",- u" • and n” A” t. • u 
mean flow equations (23) and (24) require fLthA>V a ’ 

mination of their local values Thir^> further description to provide deter- 
Urn" which Arises from sveraL'ne Tf "closure prob- 

to (18), which are believed to contain^ j-ustantaneous governing equations, (13) 

of information in the averaging process\ecLsU^ef 

resulting terms, e.g. , u",- n" • L necessitates modeling" of certain 
tion, the general problem is described fo?? ^^°\®‘inations . In this sec- 
turbulence models, variable density ef fects'^'^and ^ of two-equation 

curvature and swipl which can be corrected 


as a 


™iutL"rofthri:no^i";;:s;:Lrr""* -■> »ouna 

transport equation, Jones, 197928 ^ 




<0 

j J 


■ ‘°”j Ifr * "”i Ifi- - p (' 


i'‘”k - P 


(58) 


hav.e to bs related to lowL SrderCusuin . correlations appear, uhich 

Cartaia key assus.ptio„s regarding ^rcMve«“ /J aet. 

fuaipn term [first term on rSh?*siS, ■!»- 

brarc forms which can include* important eff.ri’ equations to alge- 

swirl, which introduce significant ani«sr.t-r streamline curvature and 

might appear in practical combustors Before^^ Reynolds stresses, as 

these so-called algebraic stress models thrstlnda^dT with 

"O'^clins approach is con^idereffirst ''“'O'" 

ap^proach are recognized at the outset, but it is deficiencies of this 

parisons for inclus ion of variable de^sitv ®°"- 

well as the above-mentioned anisotropic effects ^ t ^kk combustion, as 

have indicated man, of the difficjj^^ :Su-appj'aJh 
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Following Jones and Launder, 19721^, the Reynolds stresses are assumed 
to be linearly related to the rate of strain according to 



dtl. at! at 

■ 1*1 'Ti- -517’ * T *ij * “t 


(59) 


where E Is the mass-averaged turbulent kinetic energy, 1/2 Similarly 
the turbulent flux of each scalar quantity is related to the mean'^flow gradient 
according to 


P 


»*t 


JL 

dx. 

J 


(60) 


where is the appropriate (constant) turbulent Prandtl-Schmidt number for 
each scalar 4), which must be chosen by "optimization" with use of comparisons 
with appropriate experimental data. 


There is an implicit assumption of a unique length scale leading to the 
effective eddy viscosity From the exact transport equations for the 

Reynolds stresses, Tennekes and Lumley, 197278 have shown that this assump- 
tion IS valid if the vorticity of eddies involved in cross stream momentum 
transfer scales with the cross stream vorticity of the mean flow. Further- 
more, the shear stress has to scale with any components of the Reynolds stress 
tensor. This can occur only if there is a unique length-s.cale and time-scale 
the flow. Simple shear flows fall into this category. 

A recirculating flow such as that behind the bluff— body described in Sec- 
tion 7.1, has several characteristic length scales: the jet diameter, the 

annulus height, the thickness of the boundary layer on the bluff-body and the 
recirculation zone length. Then zero stress may not coincide with zero strain, 
Durao and Whitelaw, 1978'^. Thus the assumption of an effective eddy viscos- 
ity is questionable in recirculating flows unless the turbulence is dominated 
by mean pressure gradients. 

The isotropic eddy viscosity is given by 

p. = C ~p E^/e 

where e Is the rate of dissipation of turbulent kinetic energy. For Favre— 
averaged k and e, the following equations are solved 
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Pti 


dl 

dx. 



p'u' 


-9p _ -ry 


7®. 



'1 Pf 


3t?. p'u’ 


j-nj 


M. 


— 9x, 


(61) 



7^ __ 


(62) 


Note that in these equations, the divergence of the mean flow field 
explicity appears. From equation (22), 


90i/3xi, 


«i 

ir 

1 


P 


(63) 


which is zero only in constant (mean) density flows. 

equations have been discussed by 

Jones and Launder 197214, and many others. The constants chosen are C = 

ti£n^ i 7 / 7 ^^ ’ = 0-7. These e^ua- 

t^ons, (42) through (4?2, are, in fact, ’the variable density version through 

through the p'u"i correlations. 

h correlations pose a major modeling problem, as the above equation 

(except for the specific variable density terms) 
directly from turbulence modeling in constant density flows (including the 

above). The use of the above R-'e model in practical combus- 
tion systems, i.e., with major recirculation zones and with large dilution iet 

tro^L : pressure, r.ther then furbJleLrcon- 

trolled as in the decay of wakes, can be expected to yield realistic results 

recirculation zones and effects of high swirl 
can be refined with proposed algebraic stress models. 


i The correlations 7^7 in eouations (61) and (62) can have a significant 

iSr®;, """ Williams, 198177 , and Bilgerr 
1980 . Various models have been proposed in the literature, Jones, 197928 
but at the present the gradient flux model has been assumed in the General ’ 
Electric combustion models, viz. 
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(64) 


P'n"i =- 


9p 
cTt ax. 


This is expected to show particularly strong effects in confined flows v;here 
strong pressure fields exist, including diluti on or swirl-dominated flow 
fields. Significant work is required as these p'u"£ correlations are cou- 
pled to the mean pressure gradients in (61) and (62) and can be expected to 
have significant effects in the pressure-dominated flows in gas turbine 
combustors . 


! These terms, as well as the terms, have generally been neglected 

byi many authors, Mongia, et al., 198ll; Syed and Sturgess, 1980 ^ 8 ; and Bozsan, 
et al. , 1981^®, despite its apparent importance for combusting flows with non- 
zero pressure gradients. Furthermore, Borghi and Dutoya, 1979®^ and Schefer, 
etj al. , 198282 suggest that the conventional gradient transport expression, 

Eqi. (64), may predict the wrong sign of the turbulent transport flux, i.e., 
th at coun tergradient diffusion may occur. It is thus obvious that the role 
of j p ' u"£ in prediction of turbulent, recirculating flows needs considerable 
attention to provide more appropriate modeling and testing with detailed com- 
parisons with experimental data. Alternates to (64), including the solution 
of' transport equations for p'u"j^, have been proposed (e.g., Jones, 197928). 

The complexity of this approach requires further development and study, and 
the General Electric models presently use the simple model given by equation 
(64) . 


In addition to the above-cited density-velocity fluctuation correlations 
which couple to the mean pressure gradient, questions regarding flame generated 
tuirbulence need further study. Many investigators, e.g. Bray, 197983; Libby 
anjd Williams, 1981^^; and Bilger, 198080, have raised this issue and proposed 
alternate models for limited application (with simple diffusion or preraixed jet 
flames). For ejcample, Libby and Bray, 1978^8 multiply the eddy viscosity by a 
density ratio (p/pinitial^® ^*^8 show certain turbulence effects can be so cor- 
related at least in premixed flames. Bilger, 198080 also pointed out the pos- 
sibility of important direct effects of density on turbulence properties which 
alters the rate of generation of turbulent kinetic energy due to mean flow gra- 
dients. In an interesting study, Borghi and Dutoya, 197981 suggested that, 
due to the chemical reactions, the eddy diffusion coeffi'^ient on the fresh gas 
side of the flame should be increased, and the eddy diffusion coefficient on 
the burned gas side should be decreased significantly. This implies that the 
turbulent Schmidt number may vary significantly because of the combustion. 

I As was stressed in the previous section, the problem of inadequate turbu- 
lence modeling is further compounded by the fact that the numerical accuracy is 
not always satisfactory in a simulation process. For example, the measurements 
of| the flow in an annular jet, made by Durao and Whitelaw, 197879^ have been 
us^d for the evaluation of various turbulence models, e.g.. Pope and Whitelaw, 
19758 ^; Taylor, 1981®^; and Leschziner and Rodi, 198185^ and finite differ- 
encing schemes; e.g. , Leschziner and Rodi, 198l83; and Rastogi, 1980®°. 

A general conclusion for computational studies of separted flows is that the 
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use of finer numerical grids than most of those reported in the literature is 
necessary. The calculations done by Taylor, 1981^5^ using the standard k-e 
model and extensive grid refinements, for the turbulent flows behind bluff 
bodies show that the recirculation length is only one tenth too small compared 
with the measurements, whereas Pope and Whitelaw, 197684 reported a discrep- 
ancy of one third. Hence, care must be taken in judging the validity of a tur- 
bulence model. 


The two-equation (E and e) turbulence models such as the one described 
by equations (59) to (62) suffer other weaknesses. The assumption that the one- 
point, two-variable correlations in (59) and (60) are directly proportional to 
the mean gradients implies that the eddy viscosity is isotropic. In flows with 
significant body forces such as highly buoyant, curved, and swirled flows, 
experimental evidence shows, for example, that the ratio of turbulent Reynolds 
stress to mean flow gradient is not isotropic. Such effects have been accounted 
for in atmospheric turbulence (buoyance) problems and for flows over curved sur- 
faces (e.g., turbine blade boundary layer flow and heat transfer). In addition, 
many numerical studies have shown that the recirculation zone, either in a back- 
ward-facing step or due to high swirl, is underpredicted significantly when the 
turbulence model does not account for streamline curvature and swirl properly. 
The following discussion will be mainly concerned with the representation of 
C'^^v^ture effects on this class of turbulence closure scheme, as demonstrated 
with axisymmetric elliptic flow computations. 


Leschziner and Rodi, 198165 and Hah and Lakshminarayana, 198o87 have 
proposed algebraic Reynolds closure models for the curved turbulent flows. 
Leschziner and Rodi derived the relationship between strain and stress for the 
curved turbulent flow based on rather unrealistic assumptions such as: 



r 


au 8U. 

A • 

dn ds 


* 0 


(65) 


where Ug, U^, are streamwise (s), normal (n), and radial (r) velocity 
components. Also they include the effects of curvature on only one Reynold 
stress component as follows: 


- Ti » b> = 
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8k? k^/e^ (au /an + U /R^) n /R„] 
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( 66 ) 


and the resulting coefficient for the effective eddy viscosity 


- k^k2 


[1 + 8kJ k^/e^ OU^/an + 0^/RJ 


(67) 
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a * 1.5 

P = 0.6 


is made dependent on the local curvature through the term 



( 68 ) 


The formulation of equation (66) is inconsistent with the standard k-e model 
because higher effective eddy viscosity is predicted if kiK 2 is not arbitrar- 
ily readjusted. Furthermore, the estimation of Ug, U^, for effective eddy 
viscosity is very difficult and nearly impossible for many complicated flows; 
the formulation is restricted to two-dimensional flows with one component of 
the streamline curvature. The three-dimensional turbulent flow in a combustion 
chamber often contains swirling velocity components and the curvature has more 
than one component. Consequently many Reynolds stress components (like the 
nprmal stress along the radius of the curvature) should be handled more care- 
fully. Hah and Lakshminarayana, 198087 showed that the use of turbulence 
kinetic energy production as a primary source term in the turbulence energy 
equation [see equation (69)] makes the turbulence closure scheme insensitive 
to the streamline curvature. 


where 

P * 


U 


dt 

I 



8e ^ 


+ C 


el 



P - 




1- D-*j. 

j i*'’ 


(69) 


They proposed to use the following exact form of the source term which does 
not include mean strain rate in the energy dissipation. 
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(70) 




The non-zero Reynolds stress components are estimated using the following 
modeled Reynolds stress equation 


0 - (1 + C^) (- n^u. U - u^u. U (1 - y) - 2 0^ 

" *jl« °1 Vi " T *ij * <1 - T) “ <Y> " f »ij 


(71) 


The above model is tensor-invariant and can be used for the general three- 
dimensional flow. The model has been incorporated into General Electric ellip- 
tic flow codes and tested for a wide range of two- and three-dimensional flows 
demonstrating that the curvature effects are well represented. 


6.5 OTHER TWO-EQUATION TURBULENCE MODELS 

There are other formulations of the two-equation type of turbulence 
model. While they all retain k as one of the dependent variables, the choice 
of the second variable differs. In the k — e model it is the dissipation rate 
E|, in the k-L model it is the length-scale of turbulence L, and in the k-w model 
it is a measure of the fluctuating vorticity w. Since they all represent the 
turbulence by two scalar quantities and therefore have only one characteristic 
time and length scale, they all necessarily lead to the eddy viscosity concept 
on dimensional grounds. Because it is the very concept of eddy viscosity which 
is limiting, there is nothing fundamental to differentiate between the three 
models. On the other hand since the k-e model has been the most widely tested 
with respect to the "constants”, it is logical to recommend it over the other 
two for future use. 

Another criterion which would help to pick between the three models would 
be whether or not the modeled quantities could be measured. The turbulence 
kinetic energy k can indeed be measured (e.g., by single-point LDV) but all 
the second scalars (e, L, or w) are more difficult to measure. Typically, 
two-point LDV would be required. Thus even on experimental grounds there is 
little to choose between the models. 


6.6 TURBULENCE-CHEMISTRY INTERACTION 

In a turbulent reacting flow, the interaction between turbulence and 
chemistry is of critical importance. Three general regimes may be identified. 


147 


fool; Tf ia ™ch mailer than the kinetic time 
h!^a , 'ie ?‘<=50"«0'i" and kinetically limited. If on the other 

u ^ u instantaneously in equilibrium. Finally, If Tf ~ x„ 

both turbulent mixing and chemical kinetics must be accounted for. In this 
section several approaches to this problem are discussed leading up to the most 
appropriate. The preferred approach is the assumed shape pdf for scalar(s) 
with local moments from transport equations solved simultaneously with the 
dHcussX"'' Break-Up. ESCIMO, and Monte-Carlo models are also 

In a chemically reacting flow, the individual species mass-fractions obey 
the conservation law, equation (16) ^ 




^ (p Dy-ji) + 


(72) 


Averaging, whether conventional or density weighted, necessitates closure 

assimptions, including special treatment of the highly nonlinear chemical pro- 
auccion cemi« ^ 

d,na /°’^/®"‘o«strative purposes, the elementary reaction between species Mi 
ana W 2 at rate Kf i ^ 


**1 ^ *3 


(73) 


spLiL'^l"is*^then instantaneous production (consumption) rate (w^) of 


" ‘ V2 


(74) 


Under Favre (density) averaging, the mean rate of product 


ion IS 


- <Wj> - p kj (p Tj Tj + pT"j Y"2> 

+ P1P»£> " Y"j t j + p(pkj)'» Y"j + p(pkf)" I"jY' 


( 75 ) 


The rate constant 
simply as 


kf which is a function of temperature T may be written 
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where B is Che Arrhenius preexponenClal factor and is Che activation tem- 
perature. A direct approach to expanding this function is Co replace T by its 
mean and fluctuating components 




k- » B • T" + ( — - — ) T"* + ...1 

j2 2^4 j3 


(77) 


and Favre average Co obtain 


- B tl ...] 

2T^ P 


(78) 


where 


T 

-i << —1. 
ir* 



since typically, 


T << 




However Che higher order terms are negligible only if 

T T* ' 

-2 — << 1 

which is not true in practical flows where the temperature fluctuations are 
very large. 

In general, then 

. <W2 (T)> - Wj (<T>) [1 + F] (79) 
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where F is an expansion involving corrections between means, means and fluc^ 
tuations, and fluctuations and fluctuations. 

Hence the mean production rate is not the production rate evaluated at 
the mean temperature. Direct closure methods for the expansion series are 
generally restricted by the large activation energy 

Tjr » 1 

' ' ? 

except in near isothermal flows or in systems with small activation energies. ‘ • 

The more general closure method is based on probability density functions. 

Various models discussed below account for this "unmixedness" effect in tur- ' f 

bulent reacting flows. ! | 


Chemical Kinetics 

A chemically reacting system contains the initial reactants, the products 
and the various intermediate species produced in the course of the overall 
reaction sequence. The species conservation equations contain wj, the mass 
production rate of species i per unit volume. A general chemical kinetic 
scheme comprised of N species and R reaction steps may be written in the form 


N 

1 


\)' 






N 

5 \)' Mi j = 1. 2. ... R 
i-1 ‘J ‘ 


(80) 
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where \)'ij and\)"£j are the stoichiometric reactant and product coefficients 
for species i and reaction j. The chemical source term is given by 
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where w^j is the contribution of reaction j to the net rate 


w. . = (\)" . , 
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and kfj is the forward Arrhenius rate of reaction i k , 

depends on kfi and the equilibrium consLnhr ; ^ backward rate kfai 

dependence in the preexponential. ’ nduding a possible temperature 


B. T 


Oi ^ Ta,/T 


( 83 ) 


^=.U insl arS Inl 

times. The production rate oer nm> being in equilibrium at all 

culated. laseaad, the in.taataneoaa coilfilitwio^oftSe'tJ 
le mote efficiently solved by minimization of Gib^ Lee eJergT 

necessL^trrtt:::nr:rLe:^:af;“eJrcr(^nr 
uZni :hL“ r,;LL;jn"" 

relatively cool liner whi?LtJe cS in hrcIt^K” near the 

delay times cause a shift in the location ignition 

intermediate and pollutant species for o ^ heat releasing zones. Various 
erly accounted for with equilibri.™ k 0 and NO^. cannot be prop- 

lent flow field is very signified in'aff ^ "-^^ermore, the highly tSrbu- 
important kinetically limited species. ^ '^alue of these 


Closure at Mean Temp praf..r.m 

hinatLJL"LLL':iL„"LlnrtL' wind™''-" “'"n 

fions ( 22 ). ( 23 ), ( 61 ) aL ( 6 » fS^conL^Liorol"" 

kinetic energy and dissipation rate reSnoo^' ? momentum, turbulence 

d acnased in the aarliatLac^Ln^^LiL^ 

diffusion coefficients for each commikflM-on.i local convection- 

scaling with appropriate generalized Pranrft-/ obtained. After 

conservation equation in Newton-Raohson ro ^ numbers, the species and enthalpy 

Gaussian ellMnation. PratraL e", 9 § '“arthr/°'r" 

fractions (i = i. n) ’ , ‘ *he local species mass 

values of mass density obtained from'^tha^^^ obtained. The new field 

redetarmine rha hydroLn^L sLutilf tS“ 'hsn used to 

dynamic and thefmochemical fields is repeated hydro- 

convergence on temperature is achieved ' f°'"'"rse (at each node) 

equation is calculLed in tJL pLSLe « the'I^an 

knowingly ignores effects of fluctuation/ "'e^n temperature. The method 
onstrate finite rate kinetic effects c:;%re\^^th^^^?liS:i 

The oJ^iLralgorithf for thHhemochf ' to ideal gases, 

nonlinear source terms) and the optima/// fqoations (with strongly linked, 

tions (with weakly linked source Lrms) hydrodynamic equa- 

ing code. ^ combine to provide a rapidly converg- 


151 


By way of a demonstration, the model has been applied to an axisyrametric, 
gaseous fuel, bluff-body stabilized combustor, Roquemore, et. al., 198090. The 
hydrodynamic aspect of the flow was treated by solving the equations for con- 
servation of mass and momentum and utilizing the k-e model for the eddy viscos- 
ity. The thermochemical aspect of the flow involved 15 chemical species and 
17 reactions including NOx chemistry (Table XV). 

It must be emphasized that this demonstrative calculation employed a rela- 
tively coarse grid (30 x 27) without an assessment of numerical error. The 
objective was to demonstrate the qualitative differences between the equilib- 
rium and kinetic flow fields; also, unmixedness is explicitly ignored . A more 
rigorous simulation of this flow is presented in Section 7.1. 

The flow fields for equilibrium and kinetic flow are partially shown in 
Figures 59 through 61 and magnified in Figures 62 through 64. The flow rate of 
propane is 11.66 Kg/hr, corresponding to a jet exit velocity of about 100 m/s . 
Stagnation points from the data and the models are shown in Table XVI. 

i 

* A comparison of the location of the stagnation points from the data and 
the models shows the jet penetration to be somewhat underpredicted. Despite 
this quantitative discrepancy, which could be addressed through more exact 
turbulence and turbulence-chemistry interaction modeling, and refined numer~ 
ics, several qualitative points emerge: 

1. Both theory and data show that the effect of combustion is to increase 
the jet penetration. 

2. The flow field predicted by the equilibrium model shows stagnation points, 
while the-kinetic model (in agreement with the data which shows penetra- 
tion for m > 8Kg/hr) predicts penetration. This emphasizes the conse- 
quences of approximating the combustion chemistry with an equilibrium model 
based on the mean temperature above. Similarly the frozen flow prediction 
shows two stagnation points in agreement with the data for m ^ 14 Kg/hr. 

3. Temperature profiles immediately downstream of the jet exit. Figures 65 
dnd 66 and 0.020 m downstream. Figures 67 and 68 show that the equilibrium 
flame is narrower than the kinetic flame. The peak temperature of the 
Equilibrium flame is also higher by about 300 K. These temperature over- 
predictions result in low densities at the nodes, which then lead to 
erroneous flow field calculations. Thus, the finite rates of the major 
heat releasing reactions cannot be neglected. The kinetic temperature 
profile agree qualitatively with the data. It is interesting to note 
that the effect of random composiv^ional fluctuations on the equilibrium 
flow is also to smear out the peaks. (This is the more likely phenomena 
as shown In Section 7.1.) 

' It is repeated that this demonstrative calculation cannot be used to quan- 
titatively assess the submodels. A methodical study of this particular flow 
is reported below and the critical appraisal is then made. 
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Table XV. Kinetic Scheme for Methane Combustion. 

B. N. 3 / ^ 

= 10 J T ^ e*p <- T»7 t) / kg " Bol " » ' 
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9 OH 0 

10 H 0, 

11 OH 0! 

12 OH N 
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Figure 59. No-Fluctuation Model - Frozen Flow 
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Figure 60. No-Fluctuation Model - Equilibrium Flow. 
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Figure 61. No-Fluctuatlon Model - Kinetic Flow 
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Figure 62. No-Fluctuation Model - Frozen Flow 
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Figure 64. No-Fluctuation Model - Kinetic Flow 
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Table XVI. Stagnat io n Points from the Data and the Mode ls . 


■ jet Location of First Location of Second 

(kg/hr) Burning? Stagnation Point (ai) Stagnation Point(m) 


Data 


8 

12 

>14 

28 

No 

No 

No 

Tes 

0.070 

0.102 

Jet Penetrates 
Jet Penetrates 

0.126 

0.126 

Jet Penetrates 
Jet Penetrates 

No Fluctuation 




Model 

■ 



11.66 

No (Frozen) 

0.06 

0.083 

11.66 

Tes (Equilibrium) 

0.093 

0.133 

11.66 

Tes (Kinetic) 

Jet Penetrates 

Jet Penetrates 
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Temperature , 



Figure 6&, No-Fluctuatlon Model - Kinetic Flow. 



(ii) Assumed PDF Model 


The turbulent reacting flow model described above ignores the effects of 
local fluctuations on the flow, thus implicitly evaluating chemical source 
terms at the mean local temperature. In some situations, the kinetics are 
very rapid compared with characteristic flow times; the effect of the fluctu- 
ations then becomes more pronounced than the effect of kinetics, i.e., the 
effect of unmixedness becomes important, Hawthorne, et. al., 194991. A model 
accounting for this regime of behavior has been constructed. 

The approach is based on obtaining at each point the probability density 
function (pdf) for the minimum number of dependent variables % needed to char- 
acterize the instantaneous state of the fluid. Then the local average of any 
quantity Q is 


“ J . . . J ac"?) P (t; "x) d| (84) 


where Q (5) is thg instantaneous value of the quantity e.g. , density or tem- 
perature, and P (5; x) is the probability density function at x for the 
scalar set 5. Notjce that the mean value of Q at x is not the value of Q 
at the mean local C: 


Q(^) # Q(lf (x)) ( 85 ) 


The first approach to obtaining P (5; x) is the assumed shape model described 
next. 

Instead of writing equations for the conservation of species, the elemen- 
tal mass fractions may be treated directly. As elements are neither created 
nor destroyed in chemical reactions, their conservation equations are homo- 
geneous. With the assumption of a unity Lewis number, the element mass frac- 
tions and the; total (chemical plus sensible) enthalpy may be linked to a sin- 
gle conserved (since the conservation equatons are homogeneous) scalar, the 
mixture fraction 5 which will be sufficient to describe the state of the fluid 
in equilibrium. 5 varies between zero in one stream (the air stream) and unity 
in the other, Bilger, 197692, xhe transport equation for this conserved scalar 
is then used to derive separate transport equations for the mean (|) and the 
variance (I) , Spalding, 197193^ which are solved by the same algorithm as the 
fluid mechanical conservation equations. These equations are 
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where at = 0.7, and ut is the eddy viscosity calculated from the k-e tur- 
bulence model. Thus the first two moments of the pdf for scalar 5 are known 
at each point in the flow. 

Experimental evidence shows that the temporal fluctuations in 5 are repre- 
sentable by the 3-probability density function (3-pdf) Rhodes, 197594; 
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where a and b are given by the local moments 


.2 


b = 

I 


( 88 ) 


Realizability constraints on the fluctuations at each point x 


0 < 5"^ < I (1 - ?) 


i.e. the variance cannot be greater than that due to alternating pure air and 
fuel eddies leads to a _> 0. However, a < 1 or b < 1 cause singularities 
at 5 = 0 and 5=1 which were analytically removed before numerical convolution 
of any quantity with the pdf. 

I The assumed shape pdf model does not necessarily require equilibrium chem- 
istry. If chemical kinetics are important additional scalar(s) such as the 
rjeaction progress variable can be used to formulate a joint pdf. The number 
of moment transport equations to be solved goes up but the principle of obtain- 
ing local moments, assuming a parametric shape, deriving the local joint pdf 
and convoluting to obtain local means remains the same. 


(iii) Monte-Carlo Model for Joint Velocity-Scalars PDF 

The assumed shape pdf model is obtained from a finite number (two) of its 
moments and an assumption regarding its shape (3-pdf). Furthermore it 
describes the pdf for scalar fluctuations only so that correlations ’with veloc- 
y are not known. Thus, gradient diffusion and a separate turbulence model 
(k-e) are required. All of these restrictions can be removed by formulating 
a transport equation for the^joint pdf itself in the phase space (v. x t) 
where v is the velocity and 5 is the scalar(s) set. ’ 

This approach is under active development Pope, 198195 but is currently 
far from application to recirculating flows, primafiiy because of a step-by- 
step validation for shear flows which is now in progress partly with General 
electric collaboration. The details will not be discussed here but important 
features will be highlighted. 

The transport equation for the joint pdf P can be obtained from its defi- 
nition and the Navier-Stokes and species conservation equations. The transport 
equation has large dimensionality so that finite difference solution is intrac- 
table. Hence, the Monte-Carlo technique is used, in a form optimized for com- 
putational efficiency. In the transport equation the highly nonlinear reaction 
source terms appear in closed form. Notice that this is true also for the 
scalar(s) only pdf which is discussed next. However, in the joint pdf model, 
the convection terms both for mean and fluctuating velocities also appear in 
closed foTO. Thus there is no need to postulate gradient diffusion, which has 
recently been faulted on both experimental and theoretical grounds Schefer, 
et al., 1982 . Closure approximations are needed for the fluctuating pres- 

sure terras and the viscous and diffusive mixing. These are analogous to those 
used in Reynolds stress models. Pope, 198l95. 

described here is still a one-point pHf, i.e., it describes 
the probability of various combinations of velocities and scalars occurring 
simultaneously at a single point. Because it is a one-point , if , it carries 
no length-scale information. Multi-point pdf transport equations can be for- 
mulated but their solution would be a formidable task. Hence the one-point 
pdf model requires, for example, the e (dissipation) equation or some combi- 
nation Ck/e) to provide a length-scale. 

(iv) Monte-Carlo Model for Scalars Only 

scalar(s) pdf transport equation is sought a simpler transport 
equation results. However, the solution algorithm is still Monte-Carlo. 

As a demonstration, scalar transport in two dimensional and axisymmetric 
turbulent shear flows is calculated by means of a hybrid scheme. The aero- 
dynamic flow field IS determined from a solution of the mean mass and momentum 
equations, closed for simplicity with a one-equation model for the viscosity 
(transport equation for turbulence kinetic energy, prescribed length scale). 
Other turbulence models could easily be substituted. 
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Knowing the aerodynamic flow field, a scalar transport model is used to 
create ensembles of scalar elements in the downstream locations. The coales- 
cence/dispersion model of Curl, 196396 is then used to calculate the turbu- 
lent mixing of the ensembles. This approach does not require a turbulent 
Prandtl-Schmidt number assumption as the species (scalar) transport equations 
are not solved simultaneously with the flow equations. However, a unity Lewis 
number approximation is implicit in the reactive flow modeling. 

The raicromixing of reacting chemical species is described by the coales- 
cence and dispersion of a large number of identifiable elements comprising the 
fluid at each location in the flow. Pairs of elements are selected at random 
and assumed to: (l) coalesce (mix instantaneously at the molecular level), 

(2) remain in contact for a micromixing time, tg, and (3) disperse (separate 
into two elements of identical properties). During the micromixing time, tg, 
either equilibrium or finite rate chemistry may be assumed to occur; however, 
in view of the large number of elements at each point in the flow, special 
kinetic algorithms are needed to reduce the calculation time. Such an algo- 
rithm has been developed by Pratt, 197997 and can be used if finite rate 
chemistry is to be considered. 

The coalescence/dispersion model for scalar transport has been applied by 
General Electric to a number of mixing layer and jet experiments reported in 
the literature. Results for the free round jet experiment of Birch, et. al. , 
197898 (methane jet in air, density ratio = 0.55) and the Roshko, 197699 mix- 
ing layer (helium and neon, density ratio =7) are presented in Figures 69 
through 72. Five hundred scalar elements in each ensemble were used. The 
results are generally good with only the variance of the scalar fluctuations 
being in some error. This discrepancy has been traced to the relatively sim- 
ple turbulence models used to close the mean flow equations. 

Criticism of the coalescence-dispersion model is possible at several lev- 
els. Pope, 198195 has shown that the simple Curl model does not properly 
predict Gaussians in the limit of decaying turbulence and has suggested sophis- 
ticated biasing techniques to correct this. Perhaps more significantly, the 
Monte-Carlo scalar(s) only pdf can be questioned on the basis of what it con- 
tributes that the moment transport model does not : higher (than second 

moments of the scalar(s) pdf. Whether these higher moments are really neces- 
sary and are worth the computational expense is debatable. Furthermore, one 
of the main advantages of the joint velocity-scalar(s) pdf is lost; i.e. the 
ability to directly handle velocity correlations so that gradient diffusion 
assumptions and a separate turbulence model, e.g., k-e model, are not needed. 
For these reasons the scalar(s) only pdf approach appears unattractive. 


(v) Spalding's ESCIMQ Model 

The Eulerian approach to computing turbulent reacting flows has been cri- 
ticized by Spalding, 1978100 as being incapable of directly treating fluctua- 
tions in species concentrations and temperature. The pdf transport models are 
claimed to be restricted to simple kinetic schemes for the combustion chemis- 
try; furthermore, recent observations of coherent structures in several classes 
of flows cannot be properly described in an Eulerian reference frame. 
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Figure 69. The Normalized Decay of 
CH 4 Along the Centerline 
of a Round Jet In Still 
Air. 



Figure 70. The Mean Velocity and 

Density Profiles in the 
Self-Similar Region of 
He and N 2 With Equal 
Momentum Flux In the 
Two Streams. 



Figure 71. Decay of Scalar (CH 4 ) 
Fluctuations Along the 
Centerline of a Round 
Jet In Still Air. 


Self-Similar Region of 
He and N2 With Equal 
Mass Flux in the Two 
Streams. 
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To make up for these perceived deficiencies, Spalding, 1978i00 has pro- 
posed his Engulfment, Stretching, Coherence, Interdiffusion, Moving Observer 
(ESCIMO) model. This model combines the Eulerian and Lagrangian approaches to 
describing the flow field evolution. The model is highly intuitive and quan- 
titative and supporting evidence would be difficult to generate. However, the 
conceptual grounds for the model appear to be well founded. 

Tongues of fluid from a turbulent jet are envisioned as "engulfing" free- 
stream fluid into a jet interior. The interleaved elements of jet and free- 
stream fluid are steadily stretched by the mean shear and by turbulent fluctu- 
ations so that the thickness continuously decreases. Tpiis "sandwich" of jet 
and freestream fluid can never be separated; the material is imagined to cohere 
as it flows downstream. 

Within the sandwich, the laminar processes of diffusion and chemical reac- 
tion can occur ("interdiffusion"). These processes are analyzed in a moving 
reference frame attached to the sandwich. The processes are thus one dimen 
sional and transient and can be described by parabolic differential equations. 

Thus the ESCIMO theory consists of an Eulerian (or "demographical") aspect 
in which the "birth rate", transport, enviromment and "death rate" of the 
"'sandwiches" are described and a Lagrangian (or "biographical") aspect, which 
treats the one dimensional transient details of the "life". The results of 
the theory could then be expressed in terms of pdf's. Although this theory 
may never be detailed enough for implementation in an elliptic flow model, 
preliminary applications to simple flows such as stirred reactors is encourag- 
ing. The model has many significant physical features which can be calibrated 
(by empirical coefficients) as required.. However, there is insufficient 
direct evidence to quantitatively justify many of the key assumptions in the 
mode 1 . 


(vi) Eddy-Break Up Models 

Another semiempirical model for turbulent combustion has been postulated 
by Spalding, 197193. This so-called "Eddy Break Up" (EBU) model has applica- 
tion in situations where the local turbulence state has a dominant influence 
on the local reaction rate. The EBU model is easy to implement and has been 
used ill combustor performance codes with some success, Mongia, et al., 1979 . 

The EBU model was first formulated to predict certain features of turbu- 
lent premixed flames. The plane flame behind a rod-shaped flame holder exhi- 
bits a constant spreading angle for a variety of mixture ratios, approach 
velocities, and turbulence levels and temperatures of the inlet stream, Howe, 
1963101. Direct solution of the governing parabolic partial differential 
equations (applicable downstream of the recirculation zone) under the assump- 
tions of Arrhenius bimolecular kinetics with mean concentrations and tempera- 
tures gave flame spread rates which were of the right order of magnitude. How- 
ever, the predicted results were far more sensitive to the flow conditions than 
were the experiments. 
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To overcome this difficulty, an alternative to the Arrhenius burning rate 
was sought. The gas is imagined to be composed of discrete parcels of compo- 
sition ranging from the unburned (reactedness t = 0) to the fully burned state 
(t = 1). Then rates up to the maximum reaction rate (maximized with respect 
tot) m_ax can occur at the interface between parcels or eddies. The rate at 
which eddies are broken up — to form new interfaces — is obtained by analogy 
with the rate of decay of turbulence, 0.35 P k 3u/3y (p “ density, k * local 
turbulence kinetic energy, 3u/3y = local velocity gradient). This analogy is 
made for the equilibrium state (turbulence generation and decay rates are 
equal) . 

Two alternative rates are now available, the Arrhenius rate (niArr ^ 
the eddy break up rate (mggjj) . The EBU model caters to the sifflultaneous 
influence of each effect by actually using the hybrid rate 


^ “Art * ^ 


-1 1-1 
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The model then predicts the rod stabilized premixed flame more accurately than 
if m = niArr were used. Spalding, 1971^3 also points out various other effects, 
such as Reynolds number influence and proposes empirical formulae to account 
for them. Further developments, Spalding, 1976102, were also reported, reflec- 
ting a shift towards the ESCIMO (combined Eulerian-Lagrangian) model discussed 
However, the drawback of the EBU model is that although difficulties 
with Arrhenius kinetics at mean conditions have been identified, the proposed 
alternative approach is too empirical to implement in an overall computational 
scheme that is to be composed of "fundamental" physical submodels. In fact. 
Pope, I 977 IO 3 has demonstrated that the predictions are multiple-valued. 


Conclusion • 

In this section several modeling approaches to turbulence-chemistry inter- 
action have been discussed: 

1. Closure at mean temperature. 

2. Assumed-shape scalar(s) pdf from finite number of moments. 

3. Joint velocity-scalar(s) pdf: Monte-Carlo. 

4. Scalar(s) only pdf: Monte-Carlo and the empirical approaches. 

5. ESCIMO. 


6 . 


Eddy break up. 




It is concluded that (2) is the most tractable and accounts for the important 
physics. In the long run if computational efficiency and several detailed 
algorithms are improved the joint velocity-scalar(s) Monte-Carlo model may 
be practicable and should in principle be more accurate on physical grounds. 
The procedure outlined in (2) was used in the code GETREF and the results 
are described next. 
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7.1 MODEL ASSESSMENT AGAINST DEFINITIVE AXISYMMETRIC EXPERIMENTS FROM THE 
LITERATURE ~ 


The various submodels and techniques that are available have been ana- 
lyzed in Section 6.0 of this report. Here predictions for four axisymmetric 
experimental turbulent flows, two nonreacting and two reacting, will be pre- 
sented. The codes used are axisymmetric parabolic or elliptic (as required) 
formulations of the recommended physical models as discussed above. 

The assessment against "benchmark" data was originally intended to 
address only the three recirculating flows in Table VI. However, during 
this contract period relevant modeling and experimental work on a turbulent 
jet diffusion flame was performed (DoE Contract DE-AC04-78 ET 13146). Since 
shear flows allow evaluation of physical submodels with minimum ambiguities 
due to numerical error, and since the spectroscopic data gathered is not 
available for recirculating flows, the theory-data comparisons from the DoE 
report, Lapp, et al., 1983^, will be cited here. 

The experiments which were simulated are, in the order of the discussion 
to follow: 


1. Turbulent jet diffusion flame: Lapp, et al., 1983^ 

i , , . , 

2. L Axisymmetric isothermal slit-jet plenum/chamber : Green and 

Whitelaw, 1980 

3. Co- and counterswirled pipe flow: Vu and Gouldin, 1980®, 

4. ; Bluff-body stabilized diffusion flame: Roquemore, et al., 19809P. 


7.1.1 Turbulent Jet Diffusion Flame 

The turbulent jet diffusion flame allows assessment of an important fea- 
ture of more general flows, viz., the interaction between turbulence and chem- 
istry. Two major complicating issues, turbulence and numerical truncation 
error, are minimized. The turbulence is simpler to model because it is gen- 
erated by strong cross-stream shear; in fact, most turbulence model constants 
are derived from shear flow simulations. Furthermore, the very large radial 
gradients increase the accuracy of the gradient transport assumption. The 
issue of differencing is also simpler: being convection-dominated in the 

axial direction the governing equations are in parabolic form and axial diffu- 
sion is negligible. Thus, upwinding in the axial direction is acceptable so 
that in the radial direction grid nodes can be closely spaced. Finally, a 
forward marching solution algorithm may be used so that iteration is not 
required. 

i The experiment simulated here is a 40% CO, 30% H 2 , 30% N 2 "syngas" jet 
burning in coflowing air in a tunnel. The fuel tube is 3.2 mm in diameter 
(Figure 73). Optical access was arranged by virtue of two 0.9 m quartz win- 
dows. Data was obtained by Laser Doppler velocimetry, pulsed Raman scattering 
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Figure 73. Burning Turbulent Diffusion Jet with LDV and Pulsed Rama 
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Figure 75. Radial Profile of Axial Velocity at X/D = 
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Figure 76. Radial Profile of Axial Velocity at X/D = 100 
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Figure 80. Centerline Profiles of Jet Half-Width 
Based on Velocity (LV) and Favre - 
Averaged Mean Mixture Fraction (Raman) 
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Figure 82. Radial Profile of Turbulence Kinetic Energy at X/D = 25. 
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Versus Height Above Centerline 


Figure 84. Radial Profile of Favre-Averaged Variance in Mixture 
Fraction at X/D *25. 


measured in a cold jet into still air, Becker, et al., 1967^05; 1.3 m a 
slightly heated jet, Antonia and Bilger, 1976106 ; and 1.2 measured by Mie 
scattering in an Re » 11,200 H 2 flame, Kennedy and Kent, 198110 /. The meas- 
urements shown in Figure 84, however, do not show an Off-axis peak, which was 
expected in light of the presence of an off-axis peak in the velocity fluctua- 
tions (Figure 82). This may indicate a problem with the data rather than with 
the model and is under further study. 

Radial profiles of conventionally-averaged values of temperature are 
shown in Figures 85 and 86 at x/d«25 and 50, respectively. Compared to the 
profiles of u and % which were in excellent agreement with the experimental 
data, the modeling of temperature shows significantly larger departures. This 
is undoubtedly related to the fact that velocity and mixture fraction are 
explicit variables followed in the model while temperatures and species con- 
centrations are calculated by a further convolution with the assumed shape 
mixture fraction pdf. The mixture fraction pdf at any point in the flow is 
modeled using only two parameters u and (i"2) ^ and assuming a clipped Gaussian 

form for the turbulent part of the pdf and a correlation for the nonturbulent 
intermittent part. This is because intermittency is an important feature of 

shear flows. Since the model and experimental values of are not in 

perfect agreement, it is not surprising that values for temperature differ as 
well. 

This experiment offers the highest possible caliber of experimental data 
for a case where the numerics and turbulence aspects of the modeling are 
greatly simplified. The good agreement between the data and the prediction 
validates the model for turbulence/chemistry interactions. Further simula- 
ions including thermal and fuel-bound-nitrogen formation, Lapp et al. , 1983 , 
provide more evidence but will not be discussed here. 


7.1.2 Isothermal Axisymmetric Plenum-Chamber Experiment 

Several complicated flow phenomena are simultaneously present in the gas 
turbine combustor. To predict such a flow therefore involves models for tur- 
bulence, combustion, vaporization, radiation, etc., as well as approximations 
inherent in the numerics. The benchmark experiment of Green and Whitelaw 
i9805 addresses a "combustion-chamber-like" flow including the plenum but 
which is isothermal. Comparison with the predictions therefore provides a 
test of the turbulence model and the numerical accuracy. 

Two experiments were performed; one involved round jets while the other 
had a circumferential slit which linked the plenum with the chamber. The 
second flow is axisymmetric and is therefore computationally more convenient. 
The model was run for the axisymmetric case. 

The experimental apparatus is shown in Figure 87. Water from a constant- 
head tank enters the plenum and flows into the chamber through the slit. Laser 
Doppler Velocimetry was used to determine the mean velocity and turbulence. 
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Figure 86. 


Radial Profile of Favre-Averaged Mean Temperature at 
X/D = 50. 




The equations to be solved are those for mean continuity, axial and 
radi^il mean velocity, and turbulence kinetic energy and dissipation rate. The 
standard (isotropic eddy viscosity) k-e turbulence model was used. Boundary 
conditions were taken as symmetry about the centerline, zero normal velocity 
and logarithmic law-^of-the-wall at all solid surfaces, and zero axial gradient 
of all dependent variables far downstream. This downstream boundary was 
located approximately three inner tube diameters downstream of the gap and did 
not affect the results. The inlet boundary condition for axial velocity was 
measured and this profile was used in the calculation. Inlet conditions for 
turbulence kinetic energy and dissipation rate were obtained from assumptions 
regarding the intensity and length-scale of inlet turbulence. For example, 
the initial length— scale was taken as 3% of the annulus height. The grid was 
arrived at after successive refinements, as discussed below, and comprised 74 
(axial) X 38 (radial) nodes. 

Predicted profiles of the turbulence kinetic energy (TKE) and the mean 
axial velocity are shown (Figures 88 to 91) at four different axial locations. 
The measured data are also shown. A first inspection of these results shows a 
striking feature of this flow; the agreement on the mean axial velocity is 
significantly better than that on TKE. This implies that in the mean momentum 
equations the turbulent diffusion was not a dominant term. 

Errors in this calculation could stem from inaccurate boundary conditions, 
turbulence modeling or the discretization scheme. It is interesting to note 
that Green and IThitelaw 1980^ tried to distinguish between these sources of 
error. They studied the terms in the axial momentum equation and found that 
turbulent and numerical diffusion were comparable in the gap region. However, 
these terms were dominated by convection and mean pressure gradient. By 
refining the grid to resolve steep pressure gradients in the gap, the vortex 
formed downstream of the slit was predicted. From this, the flow was labeled 
as being "pressure dominated". 

In this study a 47 x 27 grid was first used. The vortex downstream of 
the slit was not reproduced and the profiles of turbulence kinetic energy 
showed even larger disagreement with the data particularly downstream near the 
inner wall where the experimentally observed upstream separation was important. 
Successive refinements of the grid were made without sacrificing nodes else- 
where, i.e., the later grids were nested so that prior nodes did not change. 

The vortex was then predicted. The 74 x 38 node grid used eventually is far 
more detailed than the 30 x 28 grid used by Green and Whitelaw 19805. 

The occurrence of the vortex is an important feature of this flow which 
perhaps would not be realized a priori. Therefore it is important to ade- 
quately resolve regions of steep pressure gradient in any calculation. ^ This 
requires either the time-consuming process of repeated runs or a dynamic re- 
meshing. 

Although it appears that the mean flow agreement is satisfactory, the 
reasons for the discrepancy in TKE must be addressed. These errors are pro- 
nounced in the inner chamber although far downstream the agreement improves. 
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The assumptions in the k-c model and their contribution to the total error 
have already been discussed for turbulent flows in general. In this flow, two 
particular features make the model even less accurate: the importance of the 

inner wall boundary layer and the presence of a free stagnation point. Gen- 
eration of TKE in the layer and the matching of the high Reynolds number flow 
solutions to the law-of-the-wall region are thus significant. In addition, 
the free stagnation point along the centerline represents a region where 
assumptions such as gradient diffusion and inaccurate normal-stress normal- 
strain interaction are significant. The effect of streamline curvature on the 
isotropy of stresses, and the fact that in multiple length-scale flows the use 
of an effective viscosity is not entirely justifiable also contribute to the 
total error. 

The major conclusion to be drawn, however, is that in the isothermal flow 
mean pressure gradients can dominate turbulence in determining the mean flow. 
In a combusting flow, however, the turbulence is intimately connected with the 
combustion so that any errors affect the heat release and therefore the mean 
flow. Thus, although the mean pressure gradient emerges as an importnat as- 
pect of a flow, it is not necessarily as dominating as in the present slit- 
jet/plenum/chamber flow. 


7.1.3 Isothermal Co- and Counterswirled Pipe Flow 


This sim'ilation was performed to demonstrate an improved turbulence model 
In the k-e model for turbulence, a turbulence relation between turbulent stres 
ses and rate of strain is assumed to exist. In constant density flow this 
means 




T I k + \)^ 




(89) 


The kinematic eddy viscosity Vj. which appears is obtained from the solution 
of transport equations for k and e: 




k^/e 


(90) 


where Cy is a constant (0.09) obtained from matching the model with shear 
flow data. The eddy viscosity is thus isotropic. 

In flows with significant curvature of the mean streamlines experimental 
evidence indicates that the turbulent stresses are highly anisotropic. The 
experiment of Vu and Gouldin, 1982^ provides benchmark quality data for the 


mean velocity and the turbulence field in such a flow. In the present study, 
the numerical model was run with the standard (isotropic eddy viscosity) k-e 
turbulence model and also with the algebraic stress model described in Section 
6.4. Predictions from the former model are not even qualitatively correct 
while the improved model agrees closely with the data. 

The improved model is thus very significant for gas turbine combustors 
since they are characterized by large swirl components. It may be noted that 
a previous study, Srinivasan and Mongia, 1980l0°, of this flow also identified 
Weaknesses of the standard k-e model and demonstrated an improvement based on 
a Richardson number correction. However, the current algebraic stress model 
represents a more universal approach which has now been applied to several 
flows with marked success. 

The experimental apparatus consists of two swirl generators feeding an 
inner pipe and a concentrically arranged outer pipe (Figure 92). The inner 
and annulus flows may have co-, zero, or counterswirl. The mixing region test 
section is a clear acrylic pipe. Separate blowers and plenums are used for 
the two streams. Swirl is imparted by means of vanes and is defined by the 
swirl number S 

S = exial flux of angular ■omentum 

~ [ axial flux of axial momentna ] R (91) 

where R = Ri and R = Rq are used for the Inner and outer flows, respectively. 

Measurements of the mean flow were made with a directional five-hole 
pitot probe and of the turbulence by hot-wire anemometry. Qualitative dissi- 
pation data were obtained by analog time differentiation. Further details 
regarding the experimental apparatus, the instrumentation and the experimental 
errors involved may be found in Vu and Gouldin, 1980^. This report will focus 
on the model used here and the level of agreement with the data. An important 
feature of the present study is the effort to reduce truncation errors in the 
convection term finite-difference representation so that turbulent diffusion 
is not masked by numerical error. 

Two flow cases were studied in detail, the first corresponding to co- 
swirl and the second to maximum counterswirl. In the experiment the counter- 
swirled flow produced a recirculation zone extending approximately from the 
inlet to 1.5 inner pipe diameters downstream. The occurrence of this recircu- 
lation zone provides a stringent and (in view of the gas turbine combustor 
being swirl-stabilized) necessary test of the turbulence model. 

The model was applied to this flow, first using k-e closure and then the 
algebraic Reynolds stress model described above. The computational domain, 
again chosen so as not to affect the results, extended 8.5 inner pipe radii 
4ownstream of the inlet plane. After repeated runs in which numerical diffu- 
sion errors were assessed a grid of 65 (axial) x 39 nodes was adopted for all 
the calculations. This detailed grid allows an accurate assessment of the 
turbulence model. 
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Predictions using the k-e model with isotropic eddy-viscosity are com- 
pared with the measured data for the cdswirl case (Figures 93 to 95) and the 
counter-swirl case (Figures 97 to 100). The failure to predict the co-swirl 
case is clear from the centerline axial velocity plot which shows the measured 
detay to be much slower than predicted. The effects of the swirl on the lat- 
eral turbulent diffusion of momentum are not accounted for. In the counter- 
swirl case the effect is even more dramatic since the observed center-line 
reverse flow region is simply not predicted. This is in agreement with pre- 
vious applications of the k-e model to this flow, Srinivasan and Mongia, 
1980108 . 

In situations where anisotropic eddy viscosity is clearly to blame for 
the failure of the model, several researches have developed modified two- 
equation models. These models, however, are based on somewhat empirical arg- 
uments and do not carry over to general flows. In this study a more universal 
model was used. This model has already been discribed (Section 6.4) and has 
been validated in two- and three-dimensional swirling and nonswirling flows. 

The predictions using this algebraic Reynolds stress model (Figures 101 
to 108) show that the critical recirculation zone for the counterswirling case 
is indeed predicted. The maximum reverse velocity and the downstream recovery 
rate are not exactly in agreement but the recirculation zone length is very 
close. It is repeated that this calculation is on a fine grid so that numeri- 
cal errors are minimized. 

A recent study by Sturgess, et al., 1983^®^ indicates the role of inlet 
boundary conditions in this flow, performance of the k-e and derivative two- 
equation models will certainly be affected by the boundary conditions but the 
general conclusions stand. From this simulation it is evident at least quali- 
tatively that the k-e model could be unusable in the eventual aerothermal 
model. The proposed alternative, which retains the computationally attractive 
two-equation format, has been tested without modification in a wide variety of 
flows including the present case. The success of this algebraic stress model 
suggests that it should be used henceforth. 

7.1.4 Bluff Body Stabilized Diffusion Flame 

A bluff body stabilized diffusion flame combustor has been developed at 
the Air Force Aeronautical Propulsion Laboratory (AFAPL) . This experimental 
facility was designed to provide benchmark quality data to the modeling com- 
munity. Isothermal as well as combusting flows have been studied for a range 
of fuel and airflow rates. In this section, the experiment is described, 
modeling results are presented, and the implications of the level of agreement 
obtained are discussed. Simulations were made with the code GETREF, which is 
an elliptic formulation of the Favre-averaged single-scalar assumed pdf/ 
equilibrium chemistry/k-e model. 

The bluff body stabilized combustor (Figure 109) represents an important 
step towards the "real" problem. Stabilization of the flame by recirculation 
of the hot products in the wake of the bluff body allows larger flow rates 
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Figure 94. Radial Profile of Axial Velocity at Z/R^ = 3.67 for Coswlrl Flow, k-e Model 
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Figure 96. Radial Profile of Turbulence Kinetic Energy at Z/R^ = 3.67, k-e Model 
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Figure 104. Radia l Profile of Turbulence Kinetic Energy at Z/Rj^= 3,12 for 
Counterswirl Flow, Alegebralc Stress Model. 
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Figure 105. Centerline Profile of Axial Velocity for Coswirl Flow, 
Algebraic Stress Model. 
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Another distinctive feature of the AFAPL combustor is the central fuel 
jet and its importance in determining the entire flow field. In itself the 
spread rate of round free jets requires a major modification of standard tur- 
bulence models to account for the hoop-stretching of vorticity and the conse- 
quential influence on the energy dissipation rate. Pope, 1978^1®. This 
heuristic modification cannot be implemented in a general recirculating flow. 
Thus one of the most important elements of the AFAPL flow — the jet — cannot 
be properly modeled. Of course, in the wake of the bluff body it is conceiv- 
able that the free jet physics are completely overwhelmed by other phenomena. 

The central fuel jet has another modeling difficulty. Because its diame- 
ter is small compared with the tunnel, a highly non-uniform grid in the radial 
direction is required. For reasons of numerical accuracy a small expansion 
ratio in cell-to-cell dimension should be retained. This entails a large 
number of nodes. In the axial direction an increasing cell size is used so 
that the downstream boundary condition is far enough downstream. Thus, high 
downstream cell aspect ratios result, particularly along the centerline and in 
the shear layer. On the other hand, because of the two-dimensionality of the 
mean flow, it is possible to have a large number of nodes. The gas turbine 
combustor would be characterized by a more uniform grid with fewer nodes per 
dimension. One other feature of the jet is that it remains very fuel rich up 
to the stagnation point and so may not be in chemical equilibrium. 

Finally, the instrumentation used in the AFAPL combustor can be ques- 
tioned. The LDV technique is fairly reliable if seeding levels are uniform in 
both streams, otherwise, these biases affect the predictions. This is partic- 
ularly important for predicting the location of the on-axis stagnation points. 
These locations are an important aspect of assessing the predictions. Thermo- 
couple data for the mean temperature suffers from a problem discussed earlier 
in this report: the ambiguity between Favre-averaged and conventionally- 

averaged temperature. Furthermore, temperature is an indirect measure of the 
model since the primary dependent variables are mean and variance of mixture 
fraction. Temperature is obtained after a further assumption of the pdf 
shape. Thus, the temperature comparison is ambiguous (as to which average the 
data represents) and is a simultaneous assessment of the assumption of pdf 
shape, of the assumption of equilibrium and of the accuracy of calculated 
moments. It would be highly desirable to directly compare these ingredient 
assumptions and variables directly as was done for the turbulent jet diffusion 
flame, (Section 7.1.1). 

Despite many reservations, the bluff body stabilized diffusion flame is 
an important benchmark flow on the way to reliable gas turbine combustor 
models. The AFAPL experiment is well^documented and the available and future 
data are geared towards mbdeling needs. Hence although some doubts remain, 
this flow was selected as perhaps the only such flow which could be modeled. 


Computational Domain and Boundary Conditions 

The governing equations for recirculating flow are mathematically ellip- 
tic so that downstream boundary conditions are required; these were taken to 
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?! axial gradients of all dependent variables. The 

ength of the computational domain thus becomes very important: if it is too 

short the premature tendency to attain pipe floW will affect the results. If 
It is too long then the extra nodes are essentially wa’^. ted in computing the 
flow. In this calculation the length of the domain was adopted as the short- 

affect the calculated results (approximately 5-6 

bluff body diameters). 

The number and distribution of nodes in this domain were arrived at after 
an analysis of the numerical truncation error as described below In the 
axial direction 51 nodes were used with a gradually increasing cell size. The 
lengths of successive cells increased in a geometric progression with a ratio 
1.05. In the radial direction, 36 nodes were used, clustered in the cen- 
tral jet and in the annular shear layer. This grid was held constant as a 
parametric study of the flow was made. 

The choice of inlet conditions and a parametric study of their influence 
are Included in following discussions of the AFAPL experiment. 


Comparison of Predictions With Data 

airflot^ratrof‘' 9 \“°/®^ ® •’^seline case consisting of an annular 

airflow rate of 2 kg/s and a central jet flow rate of 4 kg/hr. The predicted 
results were compared with the available data. For this stimulation ! bas^Une 
set or initial conditions was used as described below. 

. ® sensitivity analysis of the predictions to the 

initial conditions was then performed. These predictions in some cases showed 
improved agreement with the data but should not be used in any sense as the 
optimal' runs. 

Finally, with the baseline initial conditions, the central jet fuel flow 
rate was increased. Several runs were made until the predicted jet flow 

discussed above, these runs encompass 
the unsteady flow regime and so detailed studies of the influence of initial 
conditions were not performed. 


Baseline Case 

* conditions for the baseline simulation are shown in Table XVII. The 

i! turbulent pipe flow with an average exit velocity of 

about 46 m/s corresponding to the flow rate of 4 kg/hr. Turbulence kinetic 
energy in the jet exit plane is taken as 3% of the mean velocity squared while 
the dissipation rate is obtained from the inlet length scale, 3 % of the tube 

vaSaicris^!r“®" ">«ture fraction is, by definition, unity while the 




Table XVII. Flow Parameters and Initial Conditions for Baseline 
Simulation. 


Jet 


4 kg/hr 
400 K 
46.2 m/s 

3% 

3% Tube Diameter 
1 
0 


Mass Flow Rate (Air) 
Temperature 
Average Exit Velocity, 
Turbulence Kinetic Energy 
Lenght Scale of Turbulence 
Mean of Mixture Fraction 
Variance of Misture Fraction 


2 kg/s 
300 K 
47.2 m/s 

Measured Function of Radial Position 
3% of Annulus Height 

0 

0 


Mass Flow Rate (Propane) 
Temperature 

Average Exit Velocity, Uf 

Turbulence Kinetic Energy 
Length Scale of Turbulence 
Mean of Mixture Fraction 
Variance of Mixture Fraction 

Annulus 




1 ^ 

!P^ 

i m 


velocity and turbulence kinetic energy profiles are used 

e.o„ „ .r.ro^c\“-‘rttLu“e“*nLf 

from J?! the bluff body is evident 

i u ‘ LTacindis^rf „.:r?hrbi: “d, 

the second (S2) is too far. Compared with the isothermal case the 

iis 

s™Su' "f iff 

!h«« ^ comparison along the centerline (Figure 113) correctly 

shows the higher levels of turbulent kinetic energy that were 

n7T>ll"lt S"i”“on'°rSf ‘l!" “""'“““S b^psrlMents. but there it 

data as presented gives the unweighted average. 

discrepancies between the data and the predicted mean velocity and 
turbulence fields are due to several different rLsons. Some of thesfapply 
ifothermal flow as well, where also the predicted centerline reverL^ 
flow region is larger than measured, Sturgess and Syed, 19821H. 

,;. A" important feature of this flow is the central jet and the rate at 
which It spreads. This spread rate clearly affects the location of the first 

is spread rate of even a free axisymmetric jet 

8^ a notoriously difficult modeling problem. Although the~two-equSon tur- 
bulence model can correctly predict the plane shear flow, it preJic^s Jhlrthe 

trend^^iT'^h jet is higher than for the plane jet. The opposite 

trend is observed experimentally. Pope, 1978H0 has speculated thL the dif- 
ference between the two flows is the hoop stretching of vortLuy if th^ axi- 
sj^metric case. This vortex stretching coupled with the conservLion of 
angular momentum decreases the width of the vortex. Thus, the round jet has a 
greater scale reduction which leads to an enhanced spectral energy transfer 

dissipation. Modification of the k-f model tf 
process results in the correct spread rate. However^ Taylor 
198l85^has pointed out that the model cannot be employed in recitcffine 

bmhA"“L%'Jhr: rualUbhau^ bh vS!“L"fuaa- 


De vn 


as e ->0 


(92) 
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Figure 113. Centerline Turbulence Kinetic Energy 




However, the additional source term has the factor v/ e which does not go to 
zero with e. If the radial velocity is negative (which does not occur in the 
jet) this will result in nonphysical negative dissipation rates. 

Thus, although the spread rate of the jet has a significant influence on 
the entire flow, the preferred model for the free jet cannot be implemented. 

This possible source of error therefore remains in the model and can be 
removed only if the geometry of the combustor is changed. It may be noted 
that the process of vortex stretching described above is not expected to be as 
significant in the gas turbine combustor. 

Other sources of error in the turbulence model which contribute to the 
total error include: 

a„ The failure to correctly account for normal stresses and the normal 

stress-normal strain generation of turbulence kinetic energy known to 
be important at free stagnation points 

b. The assumption of isotropic eddy viscosity 

c. The implicit assumption that turbulence at any location has a 
single length-scale 

d. ; The assumption of local equilibrium of turbulence 

e. The assumption that the isothermal turbulence model closures can 
directly be used in the density-weighted model 

f. The failure to account for convection of Reynolds stresses 

g. The gradient diffusion closure for turbulent correlations.* 

These points have been discussed at greater length in the formal analysis 
of the turbulence model. 

In isothermal confined flows, the effect of inaccurate turbulence model- 
ing is often dominated by the mean pressure gradient in the momentum equations. 
However, in a combusting flow the time scale for dissipation of compositional 
fluctuations and the diffusivity of mean (|) and variance (|" 2) are obtained 
from the k and e equations. Furthermore, | and |"2 are used to generate the 
density field which "drives" all the equations. The turbulence model is 
therefore more consequential in the combusting flow. 

Contours of mean and variance (Figures 114 and 115) show that the flame 
forms a closed ellipsoidal shape as observed. The stoichiometric contour 
(? ’’ “ 0.06) gives the time-mean position of the flame which is bounded by 

those locations where approximately 
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Figure 115. Contours of Variance 


in Mixture Fraction 




Notl« .hat thia is.iuat a ph^ica, "S'a' 

actually has a continuous variation of 5, 5 sna s v 

"flame" is just an a posteriori concept. 

Siac. no data for .a.n and variance Xr 

Jhrn^'rJaUdUrrflhrg^Lent diifaalon aaau.ptlo^ and the Influence of 
pre.Lre-gradient correlations in the recirculating flow. 

a radial denaity-welghted temperature profile in f 

ahowe general agreement with themocouple Involvea 

attSmUg'lSh: :ra“pr;f'the1:?"ta"he„ here as a .^function, ^Inte»ittenc,^ie 

known to ^ause atreL!* I^the annular shear layer this 

^Uh“»^a?fwhy L prS?«ed temperature gradient is steeper than ensured. 

Til 5« ire^rt e™e“:d'Io occur anrif the Umperature fluctuations will 
L HrS' ?he overpredicti^n of centerline temperature may be due to 

assuming instantaneous chemical rat°fa;t 

Ihi:Crr;^rfppro‘rratr:kfeJ: m rich re ions 
three-dimensional sector combustor of Toral and Whitelaw, 19 

... .n;rTr“S2/S“:- 

Peclet number exceeds two. Peclet Lmbers cen- 

diagonal dominance Ite is used. Thus, numerical 

used and not that the truncation error is ^ ^circulation zone, 

'td'“aro:rsrrfec-cra^ i. in 

?^t1tn?e«ifn!dSinattd!“i;«ef:r:!'S: /agnituL of the truncation error 
must be assessed, not merely the fact that it exists. 
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4 (&z|tin^ol + Bxlcos^ol) 
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occurred in the lone region wnci.c ^Hcrrmitv to turbulent viscos- 

"?5" -iLrr'::.sre:r:L:r.uriL:.^ i. iwo»a„. a. 

any point in the flow. 

A more sophisticated analysis following McGuirk, et al., 

the various terms in the finite ditr . ; , ... .^e numerical diffusion 

truncation error. Thi» enalyeie a to I«nracy of the 

error was insignificant except in the 1 * 197957 „as applied, first 

calculation, the well-known QUICK scheme, Leonard. 1979 .was appi . 
in the momentum equations only and then in all the equations. 

The QUICK scheme has been described and analyzed earlier in this report. 
».en U%rap,Uarto the .onontun ^.-tjona there f“f 

^a:g^rir::^^o^^^:•safon:4:drrTr^orr «:ny previoua usera of the Q>,XCK 

rbi S“a :h:rrt r ricreffruar e^r^ oro:"; rr ^ . 

al., 1981^ furthermore the E equation is homogeneous and so numer- 

e equations as well. Furthermore, cne ^ enu 

to nonphysical values for denpendent variables, g , g 

Although in thla calcul.tion -lective grid deuaificatiou^euau«^^ thut^^ 
numerical diffuaion was not “ j grids, truncation error anal- 

Srr«eaterrinf"iu"he orkm^-coSaumLg aid costly. Thus, a table 
Tad ic^ta^rdtfferl^cing scheme for the convection term will be recurred. 


Variation of the Boundary Conditions 

a 1 A ovit velocity turbulence kinetic energy, and dissipation 

influence on the predictions. 
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In the baseline simulation there was no radial component to the annulus 
inlet flow. When this condition is modified by taking the radial velocity to 
be 


V = - 5% u 


(94) 


the second stagnation point moves about 6% closer to the bluff body while the 
first stagnation point is unaffected. This causes better agreement with the 
data but cannot be used to replace the baseline simulation unless radial com~ 
ponents are in fact measured. The levels of TKE are unaffected except on the 
bluff body face so that the extra strain appears insignificant. 

When the inlet turbulence length-scale in the annulus is halved to 1.5% 
of the annulus height, both the stagnation points move away from the bluff 
body. Because the inlet dissipat;- u rate is effectively doubled, the upstream 
levels of TKE are reduced, leading to less radial transport across the annular 
shear layer. By about 1.5 bluff body diameters downstream, the radial TKE 
profiles are about the same as in the baseline case. This sensitivity is 
.unfortunate because length-scale boundary conditions will entail two-point or 
autocorrelative measurements in the annulus. 

Varying the inlet turbulence level in the jet does not alter the pre- 
dicted stagnation points. The shear-generated turbulence at the edge of the 
jet essentially dominates the initial condition. 


Higher Fuel-Flow Rates ; 

On increasing the jet exit velocity the stagnation points SI and S2 were 
observed to move away from the bluff body until at a flow rate of about 8 
kg/hr the jet penetrates along the centerline. The model, with the baseline 
initial conditions, was used to study these higher fuel-flow rates. The stag- 
nation points do move outward although the underprediction of SI remains (Fig- 
ure 119). By monitoring the miuimum centerline axial velocity as a function 
of the jet exit velocity (Figure 120), the breakthrough point can be esti- 
mated. The predicted breakthrough flow-rate of about 12 kg/hr is much higher 
than the observed 7-8 kg/hr. 

The agreement between theory and data deteriorates as the jet velocity 
increases. In the transitional regime where neither the jet nor the annulus 
flow dominates, unsteady effects due to acoustic resonance, inflow perturba- 
tions or heat-release feedback phenomena may make the predictions very inaccu- 
rate. For the higher central jet velocities, the flow indeed does exhibit dis- 
tinct axial slugs of flame, named "turbules, Roquemore, et al., 1982^^-2, which 
convect downstream. The axial velocity in the turbules, which were present 
about 74% of the time, is much greatei: than in the Intervening nonluminous 
regions, Magill, et al., 1982113^ -pjjg flow is, therefore, unsteady with the 
possibility of intermittent jet penetration. 
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While the current model accounts for random fluctuations in v«>locity and 
concentrations, periodicity is not addressed. Any "coherent" structures in 
the flow thus compromise the accuracy of the predictions. Such structures can 
be acoustically excited in the annular shear layer, Corres and Pitz 19821 
and change the structure of the flow. Whether they are significant ’in practi- 
cal combustors and must therefore be included in numerical models is not 
entirely clear at present. 


Cone lus ions 

Modeling the AFAPL combustor permits certain aspects of an eventual gas 
turbine combustor model to be assessed. Although general agreement with the 
measurements is obtained, there are some significant discrepancies. These may 
be divided into two categories: some are particular to this flow field while 

the others are of general significance. 

c u in the AFAPL experiment is obtained by the recirculation 

of hot combustion products in the wake of the axisjnnmetric bluff body. Simpli- 
fication of the geometry, however, makes the entire flow very sensitive to the 
development of the annular shear layer. Periodic behavior and intermittency 
in the shear layer will affect the local temperatures and densities and the 
overall recirculation zone length. Another atypical aspect of the combustor 
is the sensitivity to the fuel jet. Proper prediction of the spread rate of 

known to require a modification to the k-e model. In the 
AFAPL combustor, an error in the jet spread affects the stagnation points and 
therefore the entire flow field. 

On the other hand, some general conclusions regarding the model may be 
arawnr 

1. The k-e model performs about as well as in isothermal flows. The 
predictions show that the peak levels of TKE occur further down- 
stream than in the isothermal flow which is supported by data from 
another experiment. 

2. Contours of mean and variance of scalar fluctuations show that the 
flame has a closed, ellipsoidal shape as observed^, Together with 
the predicted temperature field, this emphasizes the appropriateness 
of the turbulence-chemistry interaction model. 

3. Radial temperature profile predictions are perhaps affected by 
neglect of intermittency in the annular shear layer. 

4. Numerical diffusion errors have largely been avoided by concentrating 
the grid near the bluff body. In a three-dimensional flow with radial 
jets, the entire computational domain would have to be centrally- 
differenced to obtain results of equivalent quality, which does not 
currently appear practical. Thus a reliable, i.e., "wiggle-free", 
numerical scheme is desired as an alternative to first order upwind 
differencing. 
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5. The importance of having the boundary conditions in order to assess 
the performance of the numerical model is reiterated. In practical 
gas-turbine combustor flows, this involves the fuel swirl-cup exit 
flow and the dilution jet flows. 

6. Theory-data discrepancies are larger for the higher jet exit veloci- 
ties. The unsteadiness observed in the jet-dominated flow regime 
becomes increasingly significant. This effect of periodicity is not 
accounted for. 


7.2 ASSESSMENT OF THE 3-D INTERNAL FLOW MODEL 

Assessment studies of the usability of the 3-D elliptic internal flow 
model were conducted using calculations of experimental test results selected 
from the available literature, and generated as part of the overall Phase I 
program. The general purpose of these studies was to learn how to apply the 
model, assess its capabilities, and identify deficiencies in the model which 
limit its application as a combustor design tool. Detailed studies of the 
numerics and physical submodels of the kind performed with the 2-D codes were 
not done here. The majority of this effort involved calculations of dilution 
in' a crossflow. For these calculations, the data selected for comparison 
included the Walker and Kors experiments, (Reference 8), and data from experi- 
mental configurations tested as part of this program. These sets of experimen- 
tal data represented relatively simple flowfields of progressively increasing 
complexity. The remaining effort involved calculations which featured the added 
complexity of inlet swirl, fuel injection with heat release, and contour flow- 
path walls. The capability to use complex inlet swirler flows, and wall con- 
tours along the entire length pf the combustor was incorporated into the model 
as part of an ongoing IR&D program. The data selected for comparison in these 
calculations included several of the experimental configurations tested as part 
of this program. 

The initial assessment studies were conducted using data from experiment 
67 of the Walker and Kors experiments. This set of data represented a two- 
temperature trace experiment involving a singl* row of dilution jets injected 
info a crossflow at a momentum ratio of about 6. This set of data was specif- 
ically selected because the dilution hole siztt and spacing characteristics 
were similar to many conventional combustor designs. Specific details con- 
cerning these data can be found in the reference. From these data, plots of 
temperature difference ratio contours were generated for each measurement 
plane downstream of the injection point. These plotted data, shown in Figures 
121 through 124, represent averaged values over the span of two adjacent holes, 
and provide a visual presentation of the jet penetration, jet shape, and jet 
temperature dissipation. 

! 

In the first of this series of studies, an evaluation of three different 
grid networks was conducted. Each grid represented a different number of 
total nodes, but more importantly, a different number of nodes was used to 
represent the dilution hole. In each of the three grids, the (X) grid lines 
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near the aft-end boundary were compressed to more closely achieve Ax/ Ay * 1. 
This grid technique had earlier demonstrated improvement in reducing the mass 
residuals by up to 20% for the same number of interations. 

The first grid network was coarse, containing 10,395 total nodes of which 
four arranged in a square pattern were used to represent the dilution hole. 

An illustration of this grid network is presented in Figure 125. In this fig- 
ure, it appears as though the dilution hole is not properly centered. This 
results from the necessity to incorporate an extra (Z) grid line in the net- 
work to satisfy the cyclic boundary treatment featured in the code. 

The second grid network evaluated contained 21,000 total nodes of which 
25 arranged in a square pattern were used to represent the dilution hole 
(Figure 126). The only difference between this grid and the first grid lies 
in the extra detail used to represent the dilution hole. 

The third grid network evaluated contained 41,600 total nodes of which 57 
were used to represent the dilution hole. Instead of using a square arrange- 
ment, the dilution hole nodes were arranged in a more circular pattern. Due 
to the rectangular nature of the grid, a round hole shape must be approxi- 
mated. An illustration of this detailed grid is presented in Figure 127. The 
storage requirements for this particular grid approaches the storage capacity 
currently available on the IBM 3081D system at General Electric. 

Each of the threa grids were run on the Northern Research 3-D code using 
the same turbulence energy, length scale, and turbulent Schmidt number input 
values. 

Turbulence Energy Parameter (AKFAC) 3 x lO"^ 

Inlet Turbulence Length Scale Parameter (ALFAC) 2 x 10~2 

Turbulent Schmidt Number 0.9 

In the Northern Research code, the jet turbulent length scale is not an 
input, but is derived as a fraction of the effective diameter of the indi- 
vidual nodes comprising the dilution jet. In these studies this fraction was 
3 percent. These values were selected from the program sample case provided. 
The cases for the first two grids were run until reasonably converged. Rela- 
tive convergence was determined by examining changes in the solution parame- 
ters after each subsequent restart run. For the first grid (coarse), 200 
iterations were run, while for the second grid (medium), 300 iterations were 
run. In the third grid (detailed), only 200 iterations were run; it was not 
fully converged. However, the excessive expense of running the model with 
this grid led to the decision not to run additional iterations. The conver- 
gence of the solution was determined by comparing calculated results after 
successive restarts. The sum of the error on mass continuity as calculated 
v^as not used as a measure of convergence. Earlier investigations had shown it 
to be an unreliable measure of the convergence of the solution. The tablula- 
tion below shows the approximate amount of CPU (Central Processing) time and 
cost per iteration for each of these three grids. 
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Per Iteration 






CPU Time 
(Min. ) 

Costs 

($) 

Grid 

No. 

1 

(10395) 

0.27 

2.50 

Grid 

No. 

2 

(21000) 

0.64 

5.62 

Grid 

No. 

3 

(41600) 

1.28 

13.12 


Plotted results from these calculations in the form of temperature dif- 
ference ratio contours are shown in Figures 128 through 131. These figures 
are to be compared with the Walker and Kors data shown in Figures 121 through 
124. The observed lack of symmetry evident in some of the calculated results 
at the jet core represents minor imperfections in the convergence of the cal- 
culated solution coupled with the interpolating technique employed in the con- 
tour plotting routine. The jet penetration calculated with each of the grids 
is compared with that calculated from the Holdeman correlations (Reference 9) 
in Figure 132. The Holdeman correlations were derived from the Walker and 
Kors studies of dilution jet penetration and mixing. 

As observed from this figure, the less detailed grids (1 and 2), after 
200 or 300 iterations produced jet penetration levels generally in good agree- 
ment with the Walker and Kors data. The more detailed grid, (3), would 
require significantly more iterations to achieve reasonable convergence, but 
may be approaching the correlation. The shape of the jet as determined by the 
outline of the temperature difference ratio contours appears to improve with 
added grid detail. The impact of the grid of the calculated spreading of the 
jets is evident from Figure 131. The calculated jet temperature dissipation 
is considerably slower than observed from the test data. The more detailed 
the grid, the worse the problem becomes for X/D^4. Grid 2 shows more rapid 
dissipation than Grid 1 for X/D<4. Further improvement would be expected 
with Grid 3 if run to convergence. Increasing the energy level and the scale 
of the input turbulence would go in the direction of generating more rapid 
mixing, thus increasing the dissipation rate of the jet. 

In the next study conducted with the data from Walker and Kors experi- 
ment 67, the impact of the turbulence energy, turbulence length scale, and 
the Schmidt number inputs was investigated. The intent of this study was to 
determine values for these input parameters which would produce calculated 
results that closely matched the test data in terras of. jet penetration, shape, 
and jet temperature dissipation. It is recognized that searching for input 
turbulence parameters that provide agreement may be little more than finding 
a fudge factor. However, identifying the sensitivity of the results to these 
parameters may be of help in understanding the modeling problem. In any case, , 
such calculations were made and are presented in this section. To conduct 
this study, a single grid was selected. The grid selected was the medium I 

grid used in the initial study, (21,000 nodes). The values through which 
each of these three input parameters are varied are shown below. 

Turbulence Energy/Velocity^ (AKFAC) 6 x 10"^ to 1 x 10"^ 

Inlet Turbulent Length Scale/Characteristic 2 x 10“^ to 2 x 10“^ 

Dimension (ALFAK) 
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Figure 131. Calculated Results of Walker a 
Detailed Grids, X/D = 16. 
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Figure 132. Comparisons of Calculated Jet Penetration with Holdeman 
Correlation for Walker and Kors Experiment 67. 




Turbulent Schmidt Number 


0.5 to 0.9 


The jet turbulent length scale remained at about 3 percent of the node effec- 
tive diameter. 

It was determined from the results of this study that the inlet length 
scsale parameter had the greatest impact on the calculated turbulence levels, 
and the jet temperature dissipation rate. Of the range of values investi- 
gated, the best results were achieved using 2 x 10“^. Varying the turbulence 
energy parameter within the range of values investigated produced an insignif- 
icant impact on both the calculated turbulence levels and jet temperature dis- 
sipation rate. No best value was determined. Varying the turbulent Schmidt 
number within the range of values investigated produced generally minor changes 
in the calculated turbulence levels, but a moderate impact on the jet tempera- 
ture dissipation. The best results were achieved using a turbulent Schmidt 
number of 0.5. 

Using the values of these input parameters which produced the best results 
with the medium grid, calculations of the experiment 67 data were made using 
the most detailed grid, (41,600 nodes). The values selected for these three 
input parameters were: 

Turbulence Energy/Velocity^ (AKFAC) 6 x 10”^ 

Turbulence Length Scale/Characteristic 2 x 10“^ 

Dimension (ALFAC) 

Turbulent Schmidt Number 0.5 

By using the detailed grid along with the selected input parameters, it was 
intended to demonstrate that with sufficient iterations the proper jet pene- 
tration, shape and temperature dissipation could be calculated. The model, 
using the detailed grid, was run for 400 iterations consuming over 500 minutes 
processor time. The results of these calculations are compared with the 
Walker and Kors experiment 67 test data in Figures 133 through 136. 

In general, the agreement was quite good. The jet shape resembles the 
reference test data. However, the calculated results produce a shape which is 
too "round" at X/D = 2. The jet temperature dissipation for X/D > 2 matched 
well with the data. The results still showed problems near the injection 
point (X/D < 2). The calculated jet penetration is a little less than 
observed from the reference test data. Even after 400 iterations, the solu- 
tion was not entirely converged. Improvement would be anticipated with addi- 
tional iterations. Because of cost concerns, no additional iterations were 
run. The calculated results after 400 iterations predict the results from 
this test data selection reasonably well. However, the calculated results 
suggest the need to use detailed grids at the expense of large storage 
requirements and excessive run times. This would severely limit the use of 
siich 3-D codes for general use in combustor design problems. 

Following the satisfactory results obtained from Walker and Kors experi- 
ment 67 test data, 3-D calculations were performed on the Walker and Kors 
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experiment 69 test data. This data set was selected because it represented 
the same dilution pattern featured in the experiment 67 data. However, the 
jet injection momentum ratio was greater, approximately 27 compared to approx- 
imately 6 in the case of experiment 67. The higher momentum ratio is more 
representative of actual combustor conditions, particularly at the first dilu- 
tion hole or primary hole row. The teat data in the form of plots of tempera- 
ture difference ratio contours is shown in Figures 137 through 140. These 
plots represent averaged values over the span of two adjacent holes. The grid 
selected for use in the initial calculations of these data was the medium grid 
(21,000 nodes). The calculation was run for 400 iterations using the Northern 
Research code. The same turbulence and Schmidt number inputs selected from 
the studies of the experiment 67 test data were used. The jet turbulence 
length scale represented approximately 3 percent of the effective diameter of 
the 25 individual nodes comprising the jet. 

The calculated results are presented in Figures 141 through 144. Com- 
parisons between these calculated results and the measured test data reveal 
that the calculated jet penetration is considerably greater than is observed 
from the test data. The calculated jet temperature dissipation is consider- 
ably slower than observed from the test data. Adjustments to the inlet and 
jet turbulence level and inlet turbulence length scale failed to provide any 
improvement.. Results obtained from calculations of the experiment 67 test 
data had also shown insensitivity of the solution to the turbulence kinetic 
energy input parameter. One possible explanation was thought to be related 
to ^ grid dependence . 

To investigate the grid dependence, calculations of Walker and Kors 
experiment 69 test data were performed using a grid network in which the Y 
grid was refined; the X and Z grid were fixed. In this grid, the dilution 
hole was represented by nine nodes arranged in a square pattern. An illus- 
tration of this grid showing the variations in Y grid detail is presented in 
Figure 145. In the (A) version of this grid, 11 Y grid lines are featured; 
in the (B) version, 21 Y grid lines are featured; and the (C) version, 41 Y 
grid lines are featured. This overall refinement produced nearly a fourfold 
increase in the number of Y grid lines. Calculations using the (A) grid were 
run for 100 iterations. Calculations using the (B) and (C) grids were run for 
200 iterations. The extra iterations were run to achieve convergence levels 
at least as good as the (A) after 100 iterations. Results of the calcuations 
perfonned using these grids are compared with the experiment 69 test data in 
Figures 146 through 148. Surprisingly, very little change in calculated jet 
penetration is observed for the increase in Y grid. There is a slight trend 
toward improvement in the calculated jet penetration with increasing Y grid. 
One additional calculation was made using a modification to the (B) grid. 

The modification involved increasing the number of Z grid lines from 14 to 
23, while maintaining the X and Y grid. The dilution hole characteristics 
remained unchanged, nine nodes arranged in a square pattern. Figure 149 
presents an illustration of the (D) grid network. Calculations using this 
grid were intended to determine if the anticipated increase in tangential 
spreading of the jet would result in a decrease in the jet penetration. This 
calculation was also run for 200 iterations. The results are compared to 
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Figure 149. Grid 4D Used to Investigate Sensitivity of Calculated Jet 
Penetration to Grid Detail Used . 
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Figure 152. 


Comparison of Calculated Results of Walker and Kors Experiment 
from the Northern Research and Garrett Codes, X/D = 2. 
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the measured test data at all measurement planes downstream of the injection 
point. Earlier calculations had problems matching test results at measurement 
planes just downstream of the injection point. The one deficiency in the cal- 
culated results still shows a lack of early tangential spreading of the jet. 
The jet shape is too round just downstream of the injection point. The 
improvement observed in these calculations can be linked to the increased jet 
turbulent length scale, as it had been previously shown that the solution is 
insensitive to changes in the turbulent kinetic energy parameter. The results 
are consistent in that increasing the length scale reduces the turbulence 
dissipated allowing greater mixing. This reduces the jet penetration, and 
dissipates the jet core temperature more quickly. In Figures 156 and 157 
plots of the turbulence intensities and length scales generated from these 
Northern Research and Garrett calculations of the experiment 69 test data are 
compared. From this comparison it is interesting to note that despite the 
identical turbulence inputs at the inlet and the jet, the two codes generate 
significantly different turbulence levels within the flow fields. Figure 157 
|lso shows that the length scales predicted by the two programs within the 
flow field are also different from each other. However, using the modified 
Northern Research code and the medium grid, calculations of the Walker and 
Kors Experiment 67 test data were rerun. The larger jet length scale which 
had produced the satisfactory results in calculations of the Experiment 69 
test data was used in this calculation. The calculated results are compared 
to the measured test data, and to the previous calculation of this experiment. 
This comparison is presented in Figures 158 through 161. Recall that the 
previous calculation used the unmodified Northern Research code in which the 
jet length scale was internally generated as approximately 3 percent of the 
effective diameter of the individual nodes comprising the jet. It is observed 
from this comparison that the rerun calculation using the larger jet length 
scale resulted in some reduction in penetration. While the results at X/D = 
2.0 are reasonably close in penetration and temperature mixing to the measure- 
ments, at the further downstream locations, the calculated penetration is less 
than the measurements. The previous calculation using the small jet length 
scale IS in good agreement with the jet penetration observed from the measured 
test data at all locations. It is also observed from this comparison that the 
larger jet length scale more quickly dissipates the jet core temperature. 
Calculated temperature levels are in agreement with the test data at X/D > 2.0. 
The previous calculation shows considerable error at this point. However” 
further downstream the rerun calculation shows the jet almost completely mixed 
out by X/D = 8.0, whereas the test data and the previous calcullations still 
show evidence of a jet core. These results provide evidence of too much mix- 
ing. Clearly, the results of the rerun calculations using the larger jet 
length scale are not as good as those from the previous calculations. 

The inability of the calculation to produce results in good agreement 
WjLth the experiment 67 and experiment 69 data sets while using a common set of 
turbulence inputs illustrates another deficiency in these codes. Since both 
experiments represented the same hole pattern with different injection veloci- 
ties, one would expect the jet turbulent kinetic energy and length scale to 
adjust with the velocity increase in a manner maintaining relative turbulence 
intensity, and nondimensional length scale. In other words, the turbulence 
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Figure 157, Calculated Turbulence Length Scales for Walker/Kors Experiment 69 AKFAC = 6E-3 
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Calculated Results of Walker and Kors Experiment 67 Using Medium Grid with 
Larger Jet Turbulence Length Scale, X/D = 16, 
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energy and length scale inputs which produced satisfactory results for the 
calculations of Experiment 69 test data should also have produced satisfactory 
results for the Experiment 67 test data. This situation is awkward in that a 
combustor designer must have beforehand knowledge of what turbulence inputs to 
use in these codes in order to produce accurate results for the flow cases to 
be modified. 

In addition to the calculation.^ of the Walker and Kors test data, 3-D 
elliptic internal flow calculations of some of the experimental test configu- 
rations tested under this program were initiated. These calculations involved 
seven experimental test configurations (4, 5, 6, 7, 8, 13, and 17C), subject- 
ing the model to progressively more complex internal flow fields. One of these 
(17C) involved fueled injection with heat release. For details into the exper- 
imental test configurations, refer to Volume II. These calculations provided 
assessment of the capabilities of the 3-D codes with more complex experiments. 

A coarse grid network was selected to use in most of these calculations because 
of the excessive run times and costs associated with more detailed grids. This 
relatively coarse grid was selected after the study illustrated by Figures 146 
through 152 indicated that the calculated penetration of a dilution jet in 
crossflow was not strongly grid dependent. A detailed grid was developed for 
use- in modeling the test configuration with fuel injection and heat release. 
Because of the reduced storage requirements in the Garrett 3-D code, it was 
selected to run calculations of the experimental test configurations. Calcu- 
lations of the Walker and Kors data had demonstrated the close similarities 
in the results obtained from both the Garrett and Northern Research codes. 

The coarse grid selected contained 13,365 total nodes. An illustration 
of this grid network is presented in Figure 162. The grid represented a sin- 
gle cup (18") section of the 90" parallel wall annular experimental test rig. 
With this grid, multiple dilution holes could be modeled in two directly 
opposed rows on the outer and inner wall boundaries. This permitted modeling 
the two holes per cup primary row arrangement, and the three holes per cup 
secondary rqw arrangement featured in the experimental test rig. All dilu- 
tion holes were modeled using nine nodes each, arranged in a square pattern. 
Because of the lack of detail in this grid, especially in the vicinity of the 
dilution holes, the geometric area of all dilution holes on the outer wall 
boundary are the same, as are all dilution holes on the inner wall boundary. 
Therefore, different hole sizes in the dilution rows featured in some of the 
experimental test configurations are not geometrically accommodated. Because 
of the annular geometry, the contraction effect on the inner wall boundary 
grid, and the resulting impact on the geometrical sizes of the inner wall 
dilution holes is not' accommodated with this coarse grid. Directly opposed 
dilution holes intended to have the same area are larger on the outer wall 
boundary than on the inner wall boundary by an amount equal to the annulus 
radius ratio. Within the codes, checks on hole flow continuity result in 
adjustments to the injection density for holes suffering from these problems. 
The total injection momentum of the dilution jet is thus conserved. 

The turbulence energy, length scale, and turbulent Schmidt number input 
values used in all 3-D calculations of the experimental test configurations 
are listed on the following page. 
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Inlet turbulent energy parameter (AKFAC) - 6 x 10“3 
Inlet turbulent length scale parameter (ALFAC) ■ 2 x 10"^ 

Jet turbulent kinetic energy intensity » 10% 

Jet turbulent length scale * 5 x 10”^ ft 
Turbulent Schmidt Number — 0.5 

Thi, generally run for 300 iterations 

ILLsive costs ^ reasonably converged solution without 

ITr-llft ^ calculated results were transformed into plots of tem- 

contours for direct comparison with the plotted 
results from the experimental test data. 

i-hp results for experimental test Configuration 8 are compared to 

Figures 163 through 166. This configuration represented the 

DattLrf«r configurations modeled. The dilution 

pattern featured a single outer wall row of jets penetrating into a uniform 
cross flow at a momentum ratio of 50. The holes in this row featured alter- 
nating hole sizes simulating the GE/FlOl primary dilution pattern. From 
comparison with test data, it is observed that the calculation predicts a 
severe overpenetration of the larger jet. This problem is particularly obvi- 
ous in Figure 164 where the calculation predicts the larger jet to impinge on 
the opposite wall. This is reminiscent of earlier calculated results^o/the 
Walker and Kors experiment 69 test data (momentum ratio of 26) using very 

inputs. However, the turbulence inputs^sed 
JLt JLr Thf diffP P*^°<?^ced a satisfactory result for the Walked and Kors 

difference in grid density may be a contributor, but further 

efforts panned 

r the comparison that calculated results 

close to the jet injection point show considerable discrepance in iet temoera- 
ure evels. However, further downstream temperature levels show much better 
agreement With the exception of the overpenLration of the largerjet thf 
calculated result at the aft measurement plane shows many similarities with 
the measured test data. Of particular note is the calculated depressLrff 
e 0.10 contour line at 0 also seen in the measured test data. The solution 
obLinL^" additional 1,000 iterations , to see if any improvement could be ■ 

results (not shown) showed almost no change 
over results obtained after 300 iterations. cnange 

Calculated results for experimental test Configurations 4 and 5 are com- 

fered onW in^the These configurations dif- 

tered only in the jet injection momentum ratio. The dilution pattern in thes^* 

uniforrr^''”%f®^‘'“’'®‘* opposed rows of jets penetrating into a 

figuratiorr Jo^^ momentum ratio of 10 in Configuration 4 and 50 in Con- 
f g ration 5. Both rows of holes featured the same pattern of alternating 
hole sizes simulating the GE/FlOl combustor primary dilution pattern. OnL 
again evidence of overpenetration from the calculation is seen in both Lts 
of comparisons. The degree of overpenetration is more severe in the higher 
lenrf •" ° ' /^at IS particularly disturbing is that similar turbu- 

ce inputs had resulted in a calculated underpenetration of the low momentum 
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Figure 167. Comparison of Measured and Calculated Results for Experimental 
Test Configuration 4. 
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Figure 170. Comparison of Measured and Calculated Results for Experimental 
Test Configuration 4. 

Contours of 


90 

^ :a 

sS 

Si 

W r°: 


9 S 

c > 

> 

3 •«» 

*< Cf? 



IICiNAL PAQI Itt 
POOR QUALITY 










lAk 'rA*ft^3^^iM 


-''<¥•'•»: *i:s'r 



COMFItUMTIW 5 OPMSED MIHMV DILUIIOII INTO UNIFOOH CIOSSFLOU N«*5I 
1,7> INCHES DOUNSIRENH OF PNINDRV JET ROU 


fi.E.Eir. TEST C0HFTGURAT10II S AFTER 3M ITERATIONS 
MFAC=EE-3;ALFAC=2E-T:$CHIIHIT ND.=A.S:E«INJ=9ia:B*l iJ= 
GARRETT 3-D CODE 

1.7S= INCHES 0/S OF PRIOARV ROM INJECTION POINT 
■«=33.V=I5,Z=27 


vv/i /THs / 


xr/ 




v-„y 


x'\ 






i>i ■ ^^-r-r-T-r • 


* , 6 - 1 * 


PHI HVC 0.19044 


Test Data 







«*♦♦♦♦* 

-4 -2 • 2 4 6 8 

Calculation 


Figure 172. Comparison of Measured and Calculated Results lor Experimental 
Test Configuration 5. 
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Figure 173. Comparison of Measured and Calculated 
Test Configuration 5. 
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ratio Walker and Kors experiment 67 test case. Some similarities between the 
calculated results and the test data can be observed. However, in general, 
the calculated results do not accurately predict the measured test results. 

Calculated results for experimental test Configurations 6 and 7 are com- 
test data in Figures 175 through 182. These two configurations 
also differ only in the jet injection momentum ratio. The dilution hole pat- 
tern featured two pairs of directly opposed rows of jets penetrating into a 
uniform cross flow at an injection momentum ratio of 10 in Configuration 6 and 
50 in Configuration 7. This overall hole pattern simulated the primary and 
secondary hole patterns of the GE/FlOl combustor. Again these compariLns 
show evidence of calculated overpenetration especially in the region located 
at 9 clockwise from top center. However, some similarities are observed. 

The calculation reasonably predicted the hot spot centrally located within the 
annulus section in test Configuration 6 and also reasonably predicted the cold 
spot centrally located within the annulus section in Configuration 7. Pre- 
dicted temperature levels in both cases are not in agreement with the data. 

As in calculations of Configurations 3, 4, and 5, this problem is much more 
acute near the jet injection point and improves considerably at the aft 
measurement plane. While the similarities between the calculations and test 
results are of note, in general, the calculated results for these configura- 
tions do not demonstrate the desired accuracy. 


The next experimental test configuration modeled featured the same dilu- 
tion hole pattern and injection momentum ratio featured in test Configuration 
7. However, the uniform inlet configuration was replaced with a conventional 
combustor dome with counterrotating swirl cups. The swirler was modeled using 
the measured velocity data previously displayed in Figure 46 (Figure 183). 


». ^ results for test Configuration 13 are compared to the measured 

through 187. From these comparisons it is observed 
that the calculated results indicate slower mixing than is observed from the 
test data. However, at the aft measurement plane, the calculation shows the 
flow field well mixed and appearing similar to the test results, although a 
slight discrepancy in the temperature • levels is observed. 


^ The final experimental test configuration modeled represented a section 
of an actual GE/FlOl combustor. For simplicity, all dome and liner cooling 
features were blocked off. Thus, only swirler and dilution airflow were pres- 
ent. Test Configuration 17C also featured the injection of JP5 fuel with heat 
release. The overall fuel-to-air ratio was 0.016. To model this more complex 
flow problem, a model that could treat contoured walls would be desirable. 


; The ability to utilize a noncylindrical combustor wall with the stairstep 
tejchnique was briefly explored. While the technique was described in the 
available user s manuals and some of the framework was in the programming it 
wa,s not fully implemented, at least not for the aft end of the combustor. 
Figures 188 and 189 show grid and calculated results with a 3-D code for an 
FlOl-combustor shape with the stairstep technique introduced. Some difficul- 
ties remain in achieving convergence that prevented the use of this technique 
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Figure 177. Comparison of Measured and Calculated Results for Experimental 
Configuration 6. 
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Figure 179. Comparison of Measured and Calculated Results for Experimental 
Test Configuration 7. 
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Figure 180. Comparison of Measured and -Calculated Results for Experimental, 
Test Configuration 7. 
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Figure 186. Comparison of Measured and Calculated Results for Experimental 
Test Configuration 13. 
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Figure 188, FlOl Contoured Wall Grid (53X 34Y 23Z) 
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Figure 189. Combustor Calculation with Uniform Inlet (No Fuel). 
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in the assessment comparisons. However, the figures illustrate some of the 
majpr limitations of this method. First, more than half of the nodes are out- 
side of the flowpath requiring storage but not providing detail in the flow- 
path. Grid concentration is introduced in regions witnin the flowpath where 
it is not necessarily needed. Achieving grid concentration next to the walls 
is awkward. Thus it is apparent that a body-fitted coordinate system would 
have significant advantages in grid optimization and elimination of unused 
nodes. It would also permit more accurate treatment near the wall. 

However, because of the convergence problem encountered with the stair- 
step boundary introduction, the configuration was actually modeled with the 
grid shown in Figure 190, featuring parallel walls and a flat inlet. A total 
of 36,708 nodes are used of which 21 nodes, arranged in a rectangular cross 
pattern, represented each dilution hole. As with the coarser grid, all holes 
in each dilution row had the same geometric area, and the contraction effect 
on the inner wall dilution holes was present. The swirler was modeled using 
the same measured velocity data used to model the swirler in test Configura- 
tion 13. The nodes representing the swirler are shown in Figure 191. The 
fuel nozzle was modeled as having a 90" cone angle, a mean droplet size of 30 
microns, and an injection velocity of 100 feet per second. The fuel injection 
point was centrally located within the swirler. Calculations with this model 
were run for 800 iterations. Total computation costs exceeded $7,000. Calcu- 
lated results were plotted as contours of the pattern factor parameter defined 
as ; 

*^local ~ ^inlet 

"^average “ ^inlet 

Comparison between the calculated results and the measured test data at the 
combustor exit plane are presented in Figure 192. It is observed from this 
comparison that the calculation accurately predicted the location of the hot 
spots (area of high pattern factor levels). This result is perhaps fortuitous 
since the model assumptions made concerning fuel injection and parallel wall 
boundaries represent gross approximations to the actual test case. 

The presence of numerical diffusion errors was of concern for all 3-D 
calculations performed. Numerical diffusion error is the second order trun- 
cation error incurred when first order upwind numerics are used. Codes which 
use the so called hybrid scheme use central differencing (second order accu- 
rate) for cells where Peclet numbers are two or less. When Peclet numbers 
are above two, central differencing becomes unstable, and the hybrid numeric 
scheme reverts to pure first order upwind differencing. Thus, numerical dif- 
fusion errror is introduced everywhere the cell Peclet numbers exceed two, and 
the calculated mixing has little relationship to the calculated turbulence 
levels. Numerical diffusion error can be reduced through the use of detailed 
grids for codes using the hybrid differencing scheme. If enough storage is 
available, grids with sufficient detail could be used which would completely 
eliminate this error. The question is, how much grid is needed to accomplish 
this? In Figure 5 a grid network is shown which was used to obtain calcula- 
tions of the Peclet numbers for a parallel wall model of a combustor with 
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inlet swirl, and fuel injection and heat release. In this grid a total of 
38,640 nodes are used requiring approximately 90% of the available storage 
capacity. These calculations were run using the Garrett code which uses the 
hybrid differencing scheme. Plots of the calculated Peclet numbers are pre~ 
sented in Figure 6 for the (X) direction while Figures 193 and 194 show cal- 
culated Peclet numbers in the (Y) direction. In these figures the shaded 
regions represent areas within the flowfield where the Peclet number is less 
than two, or where numerical diffusion error is not introduced into the calcu- 
lations. These results demonstrate that with a grid of this detail, in most 
6f the flowfield, numerical error is introduced. As seen from these figures 
the amount of error can be substantial, as evidenced by the large regions 
where the Peclet numbers are considerably greater than two. Therefore, with 
the current grid limitations imposed as a result of available storage capac- 
rty , one cannot hope to eliminate the problem of numerical diffusion error 
with the "hybrid" scheme. 

In reviewing the calculations performed as part of the assessment of 
the 3— D elliptic internal flow codes, several outstanding deficiencies were 
exposed. The following discussion addresses these deficiencies individually. 

• While axisymmetric studies described in the previous section showed 
that grid independent results with accurate numerics could be accom- 
plished, for 3-D treatment of a combustor, numerical diffusion is 
prominent in the calculations even with the maximum storage capacity 
available. While greater storage may become available, the running 
time costs would be excessive. It would seem that the most appro- 
priate method of avoiding the numerical diffusion problem is to pro- 
vide improved numerics rather than by increasing the grid density. 

• The treatment of contoured walls by the stairstep technique wastes 
grid in regions external to the flowpath, and in unneeded regions 
within the flowpath, while grid optimization near the walls is very 
awkward. The use of a body— fitted coordinate system would greatly 
improve the ability to optimize grid for a given storage or computer 
cost limitation. 

• The over penetration of dilution jets of high momentum ratios encoun- 
tered in the initial studies resulted in considerable computer analy- 
sis investigation. Calculations of the Walker and Kors data sug- 
gested that the problem was not sensitive to grid or inlet turbulence 
boundary conditions. However, it was shown that great sensitivity 
existed to changes in the jet turbulent length scale boundary condi- 
tion. Increasing the jet turbulent length scale produced significant 
increases in turbulent kinetic energy within the penetrating jet 
flowfield resulting in faster mixing, and improved agreement with 

the high momentum ratio jet data. Calculations of the experimental 
configurations tested in this program using the inlet and jet turbu- 
lence boundary conditions which had produced accurate solutions of 
the Walker and Kors data were performed. The grid selected for these 
calculations was coarse. The results from these calculations again 
showed jet over penetration. Additionally disturbing was the slight 
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under penetration observed in the calculations of experimental con~ 
figurations featuring low momentum ratio jets. If improved numerics 
become available a reevaluation of the ability to accurately calcu- 
late dilution jets in crossflow would be appropriate. Perhaps then 
an absolute, but practical, grid requirement for satisfactory accu- 
racy could be determined. 


7*3 ASSESSMENT OF AEROTHERjMAL MODULES IN GENERAL USE 

An analysis of the GE/FlOl combustor system was performed using the com- 
bustor analytical model currently in general use at General Electric. As 
recalled, the key elements of this model include the combustor inlet diffuser 
module (DIFLO), the combustor airflow distribution module (COBRA), and the 
combustor liner heat transfer module. The purpose of this investigation was 
to assess the predictive capability of the model against the measured test 
data. The FlOl combustor system was selected to model because of the large 
amount of measured pressure and liner metal temperature data available from 
component and engine testing. 

For assessing the diffuser and airflow distribution modules, component 
test data from FlOl combustor S/N 0001 were selected from the data base. 

The operating conditions selected were: 

P 3 = 16.5 atmospheres 

T 3 = 834 K 

W 3 = 30.5 kg/sec 

For the diffuser performance analysis, the geometric features of the FlOl 
inlet diffuser system along with the selected test conditions were input to 
the DIFLO module. Total pressure losses as a function of flow were generated 
for the outer stream, the dome stream, and the inner stream. The airflow dis- 
tribution analysis was performed using the model shown in Figure 195. The 
flow areas of all the combustor airflow features were known from precise mea- 
surements made on a calibrated flow stand prior to component evaluation. Tnis 
data, along with other geometrical aspects of the combustion system, plus the 
selected operating conditions were input to the COBRA model. The diffuser 
performance curves, generated from the DIFLO analysis, were used to establish 
the combustor inlet pressures. 

. Calculated results from the DIFLO and COBRA analyses are presented in 
Table XVIII and Figures 196 and 197. For comparison purposes measured test 
data and/or experience levels are included in these figures. In general, the 
agreement in calculated versus measured levels of diffuser performance, air- 
flow distribution, and combustor pressure is very good. Past experience has 
shown that even when profile effects are simulated, the diffuser performance 
levels measured in a component test rig are not identical to levels measured 
ip an engine. This appears to be related to higher turbulence levels in the 
epgine flow discharging from the compressor turbomachinery and higher Reynolds 
number. This also accounts for a part of the small observed differences 
between the calculated results versus measured levels of diffuser performance 
airflow distribution, and pressures. ’ 









Figure 196. Conparison of Measured to Calculated Pressures and Airflow Distribution 




Figure 197. Ratio of Measured to Calculated Airflow Distribution 







The data base for the heat transfer module comparison consists of several 

-“f ''''' "O* '“"bustor. Two sots oJ Jatrwerr 

obtained at simulated takeoff conditions in reduced pressure component tests 
and two sots of data ware obtained in aotnal engine Lets. The 5«rfror 
these four tests were used to make comparisons between measured and predicted 
temperatures. The predictions were made using the standard design practice 
correlations which have been utilized at the General Electric Company for a 
number years. The combustor flow and temperature distributions wLe 
obtained from COBM analysis and served as input to calculate the heat trans- 

quantities such as boundary temperatures, heat transfer coeffi- 
cients, and flame radiation* 


Table XVIII. Summary of Analysis. 


Diffuser Performance (DIFLO Analysis) 
Outer Stream 


Experience 


Flow 

29.37% W 3 


Prediffuser Pj Loss 

l.l?% 


Total Passage Pj Loss 

2 . 88 % 

1 . 5-2.1% 

Static Pressure Rise Coefficient 

0.315 


Dome Stream 

Flow 

30.55% W 3 


Total Pijt Loss 

0.43% 

0.5% 

Static Pressure Rise Coefficient 

0.905 


Inner Stream 

Flow 

40.08% W 3 


Prediffuser P-j. Loss 

0.73% 


Total Passage P'j* Loss 

1.82% 

1.5-2. 1% 

Static Pressure Rise Coefficient 

0.460 


Mass Weighted Diffuser P.j. Loss 

1.71% 



The first set of data selected was obtained as a part of a fuels effects 
program. Reference 115. Only one panel on each of the inner and outer liners 
were instrumented in that test. Calculated results from the General Electric 
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heat transfer module for the instrumented panels were made for the simulated 
takeoff condition using these data. The selected test conditions were defined 
as : 

P 3 ■ 12,4 atmospheres 
T 3 » 811 K 

Wcomb. * 1^-65 kg/sec 
f/a « 0.0289 

Comparison of the measured and predicted values are shown in Figures 198 and 
199. Both maximum (hot streak) and average temperature distributions are 
shown in the figures. The agreement between the calculated and measured 
values is good. Shown also in these two figures and subsequent figures is the 
node network used in the heat transfer program. 

A comparison for the second set of data, using all of the Indentlcal 
assumptions on heat transfer correlations. Is shown In Figures 200 and 201 
for the simulated takeoff conditions which Is defined as: 

P 3 = 16.5 atmospheres 
T 3 = 844 K 
^comb . ~ 25.4 kg/ sec 
f/a = 0.0288 

In this configuration, the aft two panels of the outer liner were impingement 
cooled while all other panels were convection cooled by the flow of air 
between the liner and the combustor casing. All of the panels in the configu- 
ration used in the previous comparison were convection cooled. As seen in 
these figures, the predicted temperatures for most of the panels are higher 
than the measured values, while as noted above, the agreement was good for the 
fuels effects data. An investigation of the data was made to determine the 
difference in the two sets of data. Two differences were found: first, the 

circumferential instrumented locations were different, and second, the method 
of installing the thermocouple to the liner surface was different between the 
two sets of data. The measured temperature variation around the circumference 
can be as much as 111 K (200* F) as shown in Figures 198 and 199 in which the 
measurements are for many circumferential locations. It is possible that 
higher temperatures existed on the liners than were recorded in this test. 

The thermocouples were installed with thermocouple sheath parallel to coolant 
flow in the fuels effects program, while in this test, the sheath was installed 
perpendicular to the coolant flow. Figure 202 compares these two installa- 
tions. Installing the thermocouple sheath perpendicular to the coolant flow 
results in a boundary layer trip causing higher convection locally near the 
thermocouple junction. This effect can be introduced into the General Electric 
heat transfer module input. Figure 203 shows the comparison between the two 
caiculsted temperature distributions for a typical panel at sea level static 
conditions showing a calculated temperature difference of 70 K (125* F) at 
the thermocouple location. 
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Figure 202. Parallel Versus Perpendicular Thermocouple Installation. 
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Figure 203. Comparison Between Measured Linear Temperature Using 

Parallel Versus Perpendicular Thermocouple Installation. 
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Thus, when the thermocouple installation feature is included in the 
General Electric heat transfer module, agreement between the calculated and 
measured data is significantly improved. However, the thermocouple still 
seems to be located in a cool region on the liner surface. 

The predictions for the overhangs agree reasonably well with the mea- 
surements on the outer liner. Figure 201. The measurements, however, are 
much greater than the predictions on the inner liner. The thermocouple 
installations were different than for the panels and involved embedding the 
thermocouples on the hot side of the overhang. Thermocouple installation 
effects for this type of installation have sometimes allowed the thermocouple 
to read hotter than the surrounding metal. Differences in the inner and outer 
flow paths could contribute to the hotter overhang temperatures on the inner 
Ixner. The flow along the outer liner is nearly straight and at 90“ to the 
dome region, while the flow along the inner liner must- turn and accelerate as 
the passage converges. This could result in a decrement in film effective- 
ness which is not accounted for in the present analysis. An additional factor 
is that the inner liner has a view of the dome region and possibly receives a 
higher flame radiation than the corresponding point on the outer liner. These 
two factors could result in higher heat loads than assumed in the calculation. 

Comparisons of the calculated and measured data for the engine test are 
shown in Figure 204. In this test, only the forward three panels and coolant 
injection slots of the outer liner were instrumented. The comparison is shown 
for the takeoff condition defined as: 

P 3 = 25.3 atmospheres 
T 3 = 799 K 

' ^'^comb ® 45.26 kg/sec 
f/a = 0.0287 

In general, the agreement is good with the exception of panel 2 in which 
the only reading thermocouple was lower than the predicted values. 

Another set of engine data was selected from the data base. In this 
test, the aft panels of both the inner and outer liners were instrumented. 

The comparisons of the calculated and measured values are shown in Figures 
205 and 206. The comparison is for the test condition defined as: 

P 3 = 25.2 atmospheres 
T 3 =800 K 
^comb ~ 88.75 pps 
f/a = 0.032 

As seen in Figure 205, the agreement for panel 5 of the outer liner is very 
good while the predicted temperatures are much higher than the measured values 
on panel 4. Their measured temperatures are very low, being only about 122 K 
(220 F) over the coolant temperature while experience has shown that this 
panel, at times, does suffer from higher temperature problems. The agreement 


337 


LEGEND 

□ HOT SIDE NODES - AVG 
O HOT SIDE NODES - HAX 
ENGINE E70021/] TEST OATA 




Panel length, cm 



/^TlUlih . I iilf/rmn 

ffUvn 


Figure 205. Comparison of Measured to Calculated Liner Metal Temperatures 
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between predicted value and the maximum measured temperature on panel 4 of the 
inner liner is good as shown in Figure 206. The predicted values for the 
other panels are much larger than the measured values. Two of the measure- 
ments were only about 111 K (200* F) over the coolant temperature and again as 
with the outer liner experience has shown that higher temperature levels exist 
on these panels. 

The analyses to date have shown that, in general, reasonable agreement 
can be achieved between measured and calculated temperatures using the current 
heat transfer module. However, in certain regions of the combustor such as 
where the flow path converges and the flow must accelerate, improvements in 
the analysis could be made with better estimates of the velocity and gas tem- 
perature near the liners. Improved predictions of velocity and gas tempera- 
ture as expected to be provided by the 3-D elliptic program could improve the 
present heat transfer analysis. It has also been shown that care must be 
taken to examine thermocouple locations and installation methods when select- 
ing measured data sets. 

The possibility of extracting information from the INTFLOW 3-D internal 
flow model for use as input to the Heat Transfer Module was examined. This 
examination included calculations from two different grids, shown in Figure 
128, and the much finer grid shown in Figure 207. 

Pertinent results from the coarse grid model are shown in Figures 208 and 
209. In Figure 208, the temperature of the outermost node is plotted versus 
the distance from the film slot. If this temperature is to be used as a film 
or adiabatic wall temperature, it should gradually increase downstream of the 
film slot as the mainstream gas mixes with the injected air. As seen in the 
figure, the temperature at the first grid point downstream of the slot is 
already near the mainstream temperature and then decreases further downstream. 
The calculated temperature behaves more closely to a hot gas side boundary 
temperature outside of the protective cooling layer. However, inspection of 
the temperature profiles across the height of the combustor shown in Figure 
209 show a difficulty with this postulated use. The radial temperature gradi- 
ents near the wall are large and change with axial and circumferential loca- 
tions, thus illustrating that even the concept of a hot side temperature just 
outside of the protective layer is faulty. A more detailed interpretation or 
utilization of the gradients seems to be required. 

Calculations were made for the fine grid model shown in Figure 207 to 
compare directly with wind tunnel measured film protection data. The calcu- 
lated temperatures using the outer node were converted to film effectiveness 
(T mainstream - T calculated/Tmainstream - T coolant) and are plotted in Fig- 
uire 210 for several M ratios. The M ratio is the ratio of film slot mass 
velocity to the mainstream mass velocity. Also shown in the figure are curves 
for wind tunnel test data obtained at two M ratios.. At first glance it appears 
that there is reasonable agreement between the calculated values and the mea- 
sured data. Further examination shows that the trends for the M ratios are 
reversed. The wind tunnel test data show that the effectiveness decreases as 
the M ratio increases, a trend which is typical of a wide range of film slot 
designs, while the reverse trend is true for the calculated values. 
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Figure 209. Calculated Temperature Profiles from 3-D Calculations 
Using Coarser Grid. 
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Figure 210. Calculated Film Effectiveness from 3-D Calculations 
Using More Detailed Grid. 


Furthermore, it is known that significantly different fiiin effectiveness 
is found for film slot geometry changes that cannot be identified by a grid 
such as that shown in Figure 207. There is a wealth of film effectiveness 
data available at General Electric and any new film slot configuration is wind 
tunnel tested before use. Wliat is needed and not available is a method to^ 
define an effective temperature and velocity just outside the film protection 
layer from the gradients such as shown in Figure 209, that can be coupled with 
definitive wind tunnel data of film effectivenes 


7.4 SUMMARIZED ASSESSMENT 

7.4.1 Diffuser Module 

Comparison of the diffuser module calculation with measured results showed 
reasonable agreement, lliis calculation provides early in the design process a 
valuable estimate of the average pressure level in the three streams that feed 
the combustor, the dome region, outer annulus, and inner annulus passages. 
However, early in the development process these pressures are measured xn spe 
cial diffuser tests. Diffuser development modifications are made early and 
frequently remain unchanged through subsequent combustor development efforts. 
Thus the basic inputs to the Flow Distribution Module during combustor devel- 
opment are available from measurements and the Diffuser Module is not a basic 
source of error. If, however, significant three-dimensional effects were 
present in these flov^s, such as wakes from fuel nozzles or structural struts, 
they may not be adequately documented by the measurements or by the computer 
programs currently used. Three-dimensional analytical treatments could pro 
vide more detailed information than available from measurements. However, 
while some older combustor designs exhibited effects within the combustor and 
in the exit gas temperature patterns that were attributable to upstream wakes, 
all current General Electric combustor/diffuser designs show no evidence of 
wakes affecting the flow through the combustor apertures. Thus, a refined 
three-dimensional treatment is not required at this time to improve the accu- 
racy of analysis for combustor durability. However, it is expected that work 
on improved diffuser analysis methods including three-dimensional treatment 
will continue at General Electric. 


7.4.2 Flow Distribution Module 

Figure 196 shows comparisons of measurements with the calculations includ- 
ing specifically the overall check of pressure drop. There is some evidence 
of Reynolds number or approach turbulence levels affecting measured results 
because engine test results at high pressure with compressor discharge flow 
turbulence levels frequently show lower overall combustor pressure loss than 
combustor component test rigs. The flow coefficient correlations for indi- 
vidual apertures were derived generally from one atmosphere laboratory scale 
data in which turbulence levels in the approach flow are lower than expected 


engine levels. In general, the effects of two-dimensional passage flow includ- 
ing approach boundary layer flows have been adequately treated in the experi- 
ments that generated the flow coefficient data. Three-dimensional elliptic 
computer programs could be used/developed to calculate flow coefficients, flow 
angles, and flow distribution within individual apertures and used as input to 
the Internal Flow Program. Figure 87 showed a two-dimensional case where the 
external flow passage and internal flow were, in fact, more than coupled - 
they were a single solution. However, the large collection of flow coefficient 
and angle data and correlations for individual apertures available at General 
Electric allow accurate calculations without the aid of detailed flow analyses. 
Whenever new aperture shapes with unknown flow characteristics are introduced 
into designs or test hardware, additional laboratory tests of the individual 
aperture are routinely conducted for inputs into the Flow Distribution Module. 
Thus, these same flow facts are available for input into the Internal Flow 
Module. At this time, more accurate methods of determining flow coefficient 
and angle, information for the Flow Distribution Module are not significantly 
limiting the accuracy of General Electric's overall Aerothermal Model. 


7.4.3 Heat Transfer Module 

Figures 205 and 206 compare calculations of the Heat Transfer Module with 
measured metal temperature data along a General Electric combustor liner. 
Within the module, heat transfer conduction within the metal structure is cal- 
culated by routines that have long been in use at General Electric with well 
established accuracy. The convective heat transfer calculations are based on 
well established correlations in general use as modified by data and specific 
flow features that might be expected to alter the correlations. Wind tunnel 
data on special local heat transfer coefficients are measured and correlated 
where appropriate and introduced into the Heat Transfer Module calculation 
method. Film slots introduce cooling air along the hot side of the combustor 
wall which helps protect the surface from the hot combustor gases. Wind tun- 
nel data on the effectiveness of cooling air to protect the surface are avail- 
able at General Electric from a large variety of cooling slot designs, and for 
substantially new slot designs new wind tunnel data is obtained. These data/ 
correlations are utilized in the Heat Transfer Module. The hot gas side 
inputs of velocity, temperature, and flame radiation are estimated from previ- 
ous experience guided by forcing the heat transfer calculations to agree with 
metal temperature measurements in adjusting these values. The selection of 
these input values is the greatest source of error in the Heat Transfer 
Module. It is planned that the Internal Flow Module will in the future be 
able to supply more accurate values and a more detailed distribution of values 
than is available at present. The accuracy would be improved if hot gas side 
velocities and temperatures were available from a sophisticated Internal Flow 
Module even if the radiation still had to be estimated from calculations 
beginning with measured metal temperatures. It is expected that this tempera- 
ture and velocity capability will become available well before reliable flame 
radiation calculations are available from the same three-dimensional calcu- 
lated internal flow field. 
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7.4.4 Internal Flow Module 


The background understanding of the Internal Flow Module, and the mathe- 
matical and physical sources of error within the module at General Electric 
came mostly from experience with the two-dimensional code GETREF. The two 
specific three-dimensional elliptic codes assessed in this effort came from 
sources outside of General Electric. In-house knowledge of the programming 
details was very limited. Some of the sources of the mathematical and physi- 
cal errors encountered were reasonably well understood because of the experi- 
ence with the two-dimensional code. The total effort was organized and struc- 
tured based on this limited knowledge, and because of the time constraints of 
the program efficient restructuring of the effort consistent with new findings 
was not possible. The most significant finding that might have resulted in 
some restructuring was the high cost of obtaining reasonably converged solu- 
tions with the detailed grids thought to be appropriate for treating combus- 
tors with reasonable accuracy. The discovery of special difficulties with 
high momentum ratio jet penetration calculations involved more exp iration 
than originally planned and resulted in less exploration of full combustor 
analyses and comparison with gas exit temperature measurements. 

In discussing the sources of error, a list of features in the order in 
which they are encountered in conducting an analysis is followed below (see 
Table XIX). In the next section the most important possible future efforts 
are discussed with respect to priority of need. 


Table XIX. Sources of Error in 3-D Elliptic Model. 

Input-Output Capabilities 
Input Information Detail 

Fuel/Swirl-Cup-Air Insertion Angles and Velocities 
Numerics 

Turbulence Model 

Compositional Fluctuation Treatment/"Unmixedness" 

Chemical Reaction/Kinetics 

Flame Radiation Model 

Convection to Wall 

Film Slot Treatment 

Boundary Shape 

Computer Storage/Cost Effects of Framework, 
Convergence and Grid Detail 

• Misimplementation of the Modeling Concepts 
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The input and output capabilities of the three-dimensional computer codes 
are somewhat inflexible. The inflexibility is the result of trying to make 
the program easy to use for the simpler combustor calculations. One example 
is the dome input provisions. Neither code permits the introduction of a 
swirler airflow around a fuel nozzle in an annular dome. This provision was, 
however, introduced at General Electric before this present program started. 
Changes in number of grid points require changes in the array descriptions 
that were simplified at General Electric. The boundary description for non- 
cylindrical walls at the aft end involving the stairstep approach while con- 
ceptually described was not implemented in the original coding. Many other 
desirable flexibility features in the input remain uncoded. The output has 
been coupled with General Electric plotting programs. 

The importance of using the correct input details of flow, angle, and 
turbulence levels has been emphasized frequently in the literature. However, 
the knowledge of these details is often very limited. In swirl cups, for 
example, it is known that substantial changes in flow direction result from 
minor changes in the aperture shape that would not change the model geometry 
inputs at all. If this flow feature is not known for input to the calcula- 
tjion, gross errors in the resulting flowfield and fuel distribution result. 

The resulting early fuel spray distribution, even if known, is difficult to 
input properly since the individual rays, as the program is presently coded, 
must all originate from a point in the same axial plane as the swirler flpw. 
The correct coupling of the spray with the swirler airflow is difficult using 
this input approach. 

, The fuel spray treatment is limited in input detail and is not completely 
coupled with the flowfield calculations. The special problems of large drops 
penetrating into low density air, which can be important for idle pressure 
levels and below, is largely eliminated at the high pressure levels where 
durability is evaluated. Even though the present fuel trajectory capability 
in the model may not permit appropriate input, it is possible to postulate an 
input technique where the fuel is input in the appropriate regions of the 
swirler air consistent with measurements, assuming fuel that follows the air 
streamlines from the input point on. By utilizing this type of treatment, the 
detailed sources of error in the present fuel injection treatment are largely 
clircumvented . 

Artificial or numerical diffusion is a major source of error. In a com- 
bjustor calculation, this effect can be so great that the calculated mixing 
becoiaes almost independent of the turbulence model prescribing the turbulence. 
Calculated results look plausible because the general flowfield characteris- 
tjics are controlled heavily by the pressure forces rather than the exact level 
ojf turbulent mixing. With the "hybrid" scheme, the accuracy of the numerics 
ijs increased by increasing the number of nodes in the calculation until cen- 
tral differencing is possible (cell Peclet number below 2). For two-dimen- 
sional flows, increasing the node density is possible and results in second 
order accurate results which allow accurate assessment of the physical model. 
However, for the three-dimensional combustors, the grid could not be refined 
to the extent needed to permit central differencing in the critical areas. 
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The number of node points was limited by the present computer storage capac- 
ity. While future storage capabilities will presumably increase, the cost 
of running such a program becomes nearly prohibitive. The appropriate solu- 
tion would seem, to be to introduce in 5 >roved numerics. Schemes that are more 
accurate than first order upwinding, perhaps second order, but that can be 
applied at high cell Peclet numbers (greater than 2) are needed. The arti- 
ficial diffusion introduced by upwind differencing would then be eliminated. 

It should be noted that since the diffusion term is approximated by (second 
order:) central differencing it is inconsistent and unnecessary to use schemes 
of higher orders for the convection term. 

Simulation of the Air Force's bluff-body combustion experiments. Section 
7.1.4, was performed by selective grid refinement until the cell Peclet num- 
bers are below 2.0 for the critical regions. Everywhere else, analysis indi- 
cated that numerical truncation error was dominated by the other terms. This 
endorses the importance of automatic error assessment in the a posteriori 
spirit of McGuirk, st al., 1981^®. 

The k-e turbulence model itself is known to be limited in its ability to 
handle certain types of flow. The constants in the model were established 
fnbm simple shear flows. Other flows are known to require other constants for 
accurate results. As mentioned in Section 7.1.4, even round (parabolic) jets 
require different constants while swirling flows require still other adjust- 
ments. Near and downstream of stagnation points, unusually high turbulence 
levels are measured that are not predicted by the model. Briefly, several 
difficulties are introduced in simplifying the Reynolds stress equations to 
the k-c model; e.g., dropping convection of shear stress and adopting an eddy 
viscosity concept. Furthermore, this eddy viscosity is generally taken to be 
isotropic (though not in the algrebraic stress model). The dissipation rate 
equation remains another major source of error; and experimental data is not 
readily available for this quantity. 

Other two-equation turbulence models have been considered. Section 6.4, 
but it is not established that any of these are in fact better than the pres- 
ent k-e model. Further refinements have been postulated that take into account 
the multiplicity of scales that are really present. Solving the full Reynolds 
stress equations leads to higher order closure assumptions and in practice can 
be excessively expensive in computer time. The errors from the numerics treat- 
ment at present overwhelm the importance of the exact turbulence model. Hence, 
the numerics need to be improved before the turbulence model refinements are 
needed. The Algebraic Stress Model (ASM), Section 7.1.3, that was found to 
work well in swirling flows, has general characteristics in its equations. 

It deserves to be included as a refinement to the present model. Other pos- 
sible modifications to the turbulence model are not established as reliable 
improvements and are, therefore, not recommended at this time. With improved 
numerics and the present k-e model, calculated flows with useful accuracy 
should be possible. Of course, future developments in turbulence modeling 
are a valuable direction for longer term research. 


A special aspect of the turbulence is the treatment of the random compo- 
sitional fluctuations called "unmixedness . ” The usual velocity, temperature, 
and fuel/air concentration measurements are time averaged values of a level 
that is actually fluctuating rapidly with time. The calculated flow fields 
also represent time average values. As discussed in Section 6.6, a probar 
bility density function, pdf, can be used to describe the variation in values 
that constitute the mean or average value. This is especially important in 
considering kinetic or thermochemical phenomena such as the heat release. If 
the time average fuel/air ratio at a point in the flow field is stoichiometric, 
the time average temperature will be significantly lower than stoichiometric 
values even with fast kinetics (equilibrium reaction). This is because there 
will be significant over and under stoichiometric mixtures constituting the 
fluctuating levels. The pdf treats this phenomena with much better funda- 
mental physics and rigor than the eddy break-up model currently in use in the 
available three-dimensional models. The pdf method has been compared with 
instantaneous measurements in a reacting jet. Section 7.1.1, and excellent 
agireement has been demonstrated. This approach has been implemented in the 
2— P code GETREF . It is planned to adopt the pdf method for treating time 
fluctuations at General Electric in future three-dimensional elliptic codes. 

The available programs use only very crude chemical kinetics. Detailed 
kinetics coupled with the flow field calculations result in excessive computer 
time, particularly when time fluctuations are simultaneously treated. How- 
ever, for durability considerations, the combustor is operating at high pres- 
sure and the chemical kinetics proceeds very rapidly. Detailed kinetics is 
needed only if emissions or low pressure phenomena are involved. Since emis- 
sions are not within the scope of this effort concentrated on durability, no 
improved treatment is needed at this time. 

Flame radiation involves radiation from the CO 2 and H 2 O that are cal- 
culated by the model and from soot for which accurate calculation methods are 
not yet established. The flow field itself is practically unaffected by the 
inclusion of flame radiation in the model. The lack of an adequate flame 
radiation model in the codes available at General Electric resulted in exclud- 
ing the evaluation of flame radiation from these present efforts. It is 
planned to postpone at General Electric efforts to provide a flame radiation 
model coupled with the three-dimensional Internal Flow Model until improve- 
ments are accomplished on other aspects of the model. 

The calculation of the convection to the walls does not provide a good 
calculation of metal temperature but the method is adequate for defining the 
^nternal flow field away from the walls. The cooling film treatment even 
yith a fine grid near the film does not produce the same wall protection char- 
acteristics that are found from measurement. The current coding does not per- 
mit the slot itself to be divided into a fine grid; the first grid point is 
intended to be at a specific fraction of a slot height outside of the slot. 

With a separate program, perhaps two-dimensional grid refinement would be pos- 
sible. However, detailed differences within the slot itself are known to 
effect the protective ability of the cooling air, but this finer grid by 
itself cannot resolve these differences. The extraction of flow conditions 
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at the boundary between the film protection layer and the hot combustor 
gases would permit the use of existing correlations of film protection char- 
acteristics. However, as discussed in Section 7.3, this boundary is complex, 
involving temperature, velocity, and turbulence gradients. It is planned to 
develop a technique at General Electric for identifying "effective" boundary 
conditions for input to the Heat Transfer Module. 

The ability to utilize nonclylindrical combustor walls with the stairstep 
technique has difficulties as discussed in Section 7.2. Grid line concentra- 
tibns can exist in regions where they don't benefit the calculation. Also, 
many node points are outside of the flow field. This causes a given number of 
nodes to produce a less accurate answer than for a combustor with cylindrical 
walls and is, therefore, a significant source of error. A body-fitted coordi- 
nate system would eliminate this problem while simultaneously permitting a 
useful concentration of grid next to the wall. 

Computing time costs as a source of error are a major consideration. 

Errors due to inadequate grid cannot be simply corrected by adding more grid. 
The storage limits of available computers are easily saturated, and detailed 
grids result in excessive computing time and costs. Looking back in history, 
the storage capacity of computers has increased at a rapid rate and computa- 
tion costs have come down. It is expected that these trends will continue. 
HoVever, the desires for increased grid are extreme. For example, an increase 
in grid by a factor of three in each direction would seem to be a moderate 
refinement, but this involves a factor of three cubed or 27 in storage require- 
ment. Furthermore, the number of iterations required to achieve a converged 
solution increases for larger node arrays. With the present costs of a con- 
verged solution for detailed node distributions utilized in this program 
exceeding $5000 per case, significant increases in nodes are not practical 
with the present codes. Efforts are needed to provide more accurate answers 
with the same number of nodes as well as to provide significantly faster algo- 
rithms or other faster techniques for achieving converged solutions. 

In addition to errors inherent in the mathematical and physical modeling 
concepts as discussed above, errors can enter the calculations due to mis- 
implementation of the concepts. This can occur because of mistakes or poor 
judgment in the actual program coding, and can also occur when the user incor- 
rectly inputs to the programs. This later possibility is increased when the 
user is not the writer of the program, perhaps even working in a different 
company, resulting in misunderstanding of parameters and input needs. More 
coinplete users-manuals, more like programmers manuals, can help. Careful 
checking of output details and features against physical expectations and 
previous calculation including those done by others and with other computer 
programs are also valuable guidance. 


Needs for Benchmark Quality Data 

The sources of error in future models will require further conq>arisons 
with benchmark quality data. Additional new experiments are appropriate. In 


addition to more experiments similar to those that have been conducted in the 
past, some experiments should give special attention to obtaining as complete 
a documentation of the flow field under study as possible. The existing data, 
including the data selected and generated in this program, frequently do not 
include as much information as the modeler generally desires for his compari- 
sons. For example, in addition to velocity, temperature, and turbulence inten- 
sity measurements in detail, the modeler who is trying to develop improved 
physical submodels would also benefit from measurement on the turbulence scale 
and the nature of the fluctuations, which are seldom obtained. Also, the 
inlet conditions to the experiment need to be correctly documented. In real 
combustors, little data exists on the internal flow details of any kind. How- 
ever, the use of models in design will require inlet flow conditions. Experi- 
ments to provide detailed information should include noninterfering techniques 
to document those features of the flow that are expected to be affected by 
intrusive techniques. In the comparisons conducted in this work, calculations 
of turbulence intensity and scale were not always available and more complete 
data on the measured turbulence would have assisted in identifying the sources 
of error. 

One specific area where additional data are critically needed in order to 
obtain reasonably accurate documentation for flow fields is in the region of 
the fuel injector and swirl cup. The flows emerging from this region are not 
properly defined at present by the model from simple representations of the 
geometry. This region can be experimentally explored to a considerable extent 
without the presence of the downstream portion of the combustor and also with- 
out reaction. This can simplify the experiment while still supplying valuable 
data on the character of the fuel and airflow in this region. 
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8.0 CONCLUDING REMARKS 


8 . 1 SUMMARY 


For comparison with the aerothermal model, data were selected from 
General Electric's data base and from the literature and new experimental 
data were obtained. The new experimental data itivolved temperature distribu- 
tion data from a series of combustor-like flows of progressively increasing 
flow complexity. Comparisons of modules within General Electric's overall 
aerothermal model were conducted. 

The comparison showed that the modules already in use are resonably 
accurate. The principal weakness in the Heat Transfer Module is in the 
accuracy of the hot gas side input. Also, the currently used methods do not 
include a module for accurate analytical estimation of the combustor exit gas 
temperature. 

The mathematical and physical principles of the internal flow model were 
investigated in one- and two-dimensional studies of various flows for which 
benchmark" quality data were available. These studies addressed the "numer- 
ical diffusion" error due to the use of first-order accurate discretization 
schemes as well as turbulence and turbulence-chemistry interaction phenomena. 

The one-dimensional prescribed uniform velocity convection-diffusion equa- 
tion was studied with and without source terms to ascertain the performance of 
various difference operators applied to the convection term; diffusivity was 
taken as a known constant. As is well known, differencing of the odd derive-^ 
tive can lead to instability if central differencing is used when the grid is 
coarse (local Peclet number less than two). One-sided upwind differencing 
provides an unconditionally stable finite difference representative of the 
convection term but is accurate only to first order in the grid spacing. This 
IS one order less accurate than the centrally-differenced diffusion term 
(which, being an even derivative does not have the stability problem) and leadk 
to artificial diffusion. Thus, the convection term differencing operator must 
compromise between formal accuracy and stability. One-dimensional prescribed 
nopunifona velocity flows and two-dimensional uniform velocity flows were also 
modelled. All these studies indicated the need for a second-order scheme 
(e.g., second-order upwinding or QUICK), rather than a first-order scheme 
(e.g., skew first-order upwinding which recognizes grid to streamline skewness 
but IS still overly diffusive). Direct central-differencing of the convection 
terra in three-dimensional flows is precluded by the excessive grid needed to 
ensure local Peclet numbers below two. 

In the assessment of the physical submodels in two-dimensional recircu- 
lating flows, however, the grid density was increased until local Peclet 
numbers allowed central differencing. Such numerical error as existed was 
quantified. With this sound numerical basis the model was applied to three 
recirculating flows. A parabolic formulation - for which the numerical issue 
IS not relevant - was also used to simulate a non-premixed jet flame for which 
high quality optical data were available. These studies indicated that the 
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suggested advanced physical models are necessary. Specifically, a probabilis- 
tic treatment for the random compositional fluctuations which lead to nonlinear 
turbulence-chemistry interaction is required. Furthermore, the isotropic eddy 
viscosity k-c model was shown to be inadequate for flows with strong swirl 
components; an algebraic simplification of the full Reynolds stress equations 
which retains the two-equation format was suggested and validated. 

In utilizing the 3-D calculations in INTFLOW difficulties were encountered 
in obtaining good comparisons with measurements of dilution jet penetration and 
mixing data. Variations of grid detail and input quantities to the calcula- 
tions were examined to determine their effect on the solution. It was found 
that the higher momentum ratio jets were more sensitive to these inputs than 
lower momentum ratio jets which are less affected in their turning by the jet 
shape development from turbulence and more by the pressure field effect on the 
early jet. Calculations done with one set of parameters, but with both of the 
3-D codes within INTFLOW, showed that similar penetration and mixing were cal- 
culated, even though the local calculated turbulence levels were different. 

The comparison of calculations for the dilution jet experiments conducted at 
General Electric are close in character to the measureriients but not good enough 
to be a useful design tool at this time. 

This effort did not establish the grid density needed to adequately treat 
dilution jets. To stay within reasonable computing time and storkge demands, 
it is not possible to calculate combustor dilution patterns accurately with 
the available codes. Tlie grid requirements for future 3-D codes with improved 
numerics could not be determined in this work; studies conducted with a 3-D 
code haying the desired numerics improvements will be required to make this 
determination. 

Comparison of calculations with an experiment that included fuel injec- 
tion and heat release and a combustor dilution pattern, showed that the calcu- 
lated gross temperature pattern in the gas of the combustor exit had consider- 
able similarity to the measured values. Further analytical and experimental 
study is still needed, however, to define an adequate method for defining the 
input for the initial fuel and air distribution. 

The computer running time costs of obtaining reasonably converged solu- 
tions was very high; thus future modeling improvements should strive for 
faster convergence techniques to increase the size of the maximum affordable 
grid along with alternate numeric schemes to provide greater accuracy for 
whatever grid is adopted. For inputting the dome air and fuel injection 
details, more experimental knowledge of the early flow development is needed. 

8 . 2 PRIORITIZED RECOMMENDATIONS 

Table XX presents a prioritized list of the main recommendations for 
treating the sources of error in the current 3-D elliptic internal flow code 
INTFLOW. Tlie list includes two paths: needed basic capability and applica- 

tion to durability. The efforts in the second category should begin whenever 
it is decided that the needed basic capability is progressing fast enough to 
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justify work on features to permit the application to the design/development 
guidance process. These recommendations are based on a methodical study of 
2~D flows and model problems and also on a more pragmatic application of the 
available 3-0 code. 

fh® first substantial improvement in capability can come from revisions 
in the basic code framework. A body-fitted coordinate system would permit 
needed improvemehts in the accurate treatment of the boundary, allow desired 
grid concentration near the walls, permit more flexible treatment of film slot 
heights, eliminate grid points outside of the flow field with their computer 
storage demands, and increase accuracy for all varying grid sizes throughout 
the flow field. Methods of accomplishing a body— fitted coordinate system can 
include a basic nonothogonal coordinate framework without any transformation, 
or a basic coordinate transformation technique. The latter seems most attrac- 
tive at present. 

Improved numerics shoiild be incorporated, particularly a technique that 
will improve first order upwind differencing by extending regimes in which 
second-order accurate schemes can be applied. Several methods have been dis- 
cussed, but further study and evaluation would be recommended before select- 
ing the most appropriate technique. 

Since running costs are a serious limitation, improved convergence speed 
is also needed. Convergence speed is affected somewhat’ by the numerics adopted, 
but solution algorithms should be adopted to give reduced running costs and, 
when larger computer storage becomes available, to permit increased grids. 

The present pressure-correction equation approach to handling continuity is 
another serious drawback in that it requires a staggered grid (with attendant 
increase in storage requirements) and also can slow down the convergence rate. 

The second item in the needed basic capability list is improved turbu- 
lence and turbulence/chemistry interaction treatment. This may be of less 
overall impact than the first item, but at least the first two improvements 
recommended in this second item have already been implemented in 2-D at General 
Electric. The algebraic stress model for dealing with the turbulence anisot- 
ropy in swirling flows has been demonstrated and deserves to be included as an 
algebraic correction to the turbulence model. 

The pdf, probability density function, treatment of random compositional 
fluctuations has been thoroughly evaluated. It is far superior in the charac- 
ter of its fundamental physics to the eddy break-up model in the current 3-D 
programs. This change is a very worthwhile improvement. 

An appropriate change to the basic two-equation turbulence model to 
handle all flow regimes has not been identified. Reynolds stress transport 
equations and improved dissipation rate equations may be the only option. 
Timewise, such a solution may be a long way into the future. 

The third item calls for improvement in the boundary flow input treat- 
ments. Inputs based on the physical geom.etry are not adequate. The flows 
leave the entry regions on angles and with turbulence levels that are affected 
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by small geometry features that are not directly treated by the models. Accu- 
rate inputs for these features must come from; empirical data. An extensive 
data base exists at General Electric, but additional data will ultimately be 
needed, and methods of translating this empirical information into inputs to 
the model need to be developed. Detailed analysis with a 3-D model of the 
flow through a specific aperture has been considered as a technique for gener- 
ating the inputs for the larger Solution, but the basic need would still 
involve a significant new data base. ^ 

The last item in the list of needed basic capabilities is an improved 
chemical kinetics treatment. This will become particularly important if the 
sme model is to be used for emissions or blowout calculations. For dura- 
bility alone its impact is small and, therefore, has been placed far down the 
list. Expanded joint pdf methods have, however, been given considerable atten- 
tion at General Electric, development concepts are available. 

On the application to durability the two durability design needs of liner 
heat transfer inputs and pattern factor calculations are listed. For both, 
a method for fuel and swirler air insertion is needed. As in the third item 
above, this involves accumulation of new data as well as an accompanying ana- 
lytical study on how to use the data. 

The boundary between the film cooling protection and the hot gas side 
flows is not distinct. Temperature, velocity, and turbulence gradients are 
all present. However, similar techniques that would be used to introduce 
finer grid near the wall to calculate through the film can also be used to 
define effective" hot gas side conditions at the film protection boundary. 

This will permit the use of existing wind tunnel data on the effect of film 
protection slots of different designs. This technique needs to be defined 
and programmed. 

Methods of treating flame radiation can be developed at a later time. 

New data is being accumulated on flame radiation in ongoing programs and this 
data could be used in calibrating any future calculation methods. The final 
flame radiation calculations are expected to be a hybrid calculation not inte- 
grated with the INTFLOW iterations but utilizing the final flow field. 

The exit gas temperature pattern factor design tool is conceived as a 
perturbation calculation at this time. While the absolute temperature pattern 
may not be accurate, much greater accuracy would be expected for the perturba- 
ation effect due to a change in dilution hole pattern. This tool would be 
developed initially along this principle. 
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APPENDIX 

NOMENCLATURE 


a, b 

- Exponents 

in B depending on ^ and 

Ae 

- Flow area 


B 

- Arrhenius 

preexponential factor 

Cl> C 2 , Cy, 

- Constants 

in turbulence model 


S 

d 

D 


f 

y 


g 

8i 

Gw 

h 


hi 

hQ 

I 


~ Specific heat at constant pressure 

- Drop diameter 

- Mass species diffusion coefficient 

- Diffuser efficiency 

- Friction factor 

- Ratio of specific heats 
” Gravitational constant 

■■ Gravitational acceleration vector 

- Jacobi iterative matrix 

- SOR iterative matrix 

- Specific mixture enthalpy 

- Specific enthalpy of species i (including heat of 

formation) 

- Heat of formation of species i 


I 


k 


hf(b) 



- Identity matrix 

- Wave number 

- Turbulence kinetic energy 

~ Forward (backward) chemical rate constant 

- Forward Arrhenius rate of reaction j 
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I''-'-'-- 


L 

L/D 

Mi 

M 

N 

Nu 

P 

HK, h 
Pe 

^eAx 

Pr 

Ps 

Pt 

Q 

R 

Rc 

Re 

Sc 

s» n, r 

t 

T 

Tc 

Td 

Ts 


- Lower triangular matrix 

- Length to hydraulic diameter ratio 

- Molecular weight of apecies i 
■“ Number of species 

“ Number of nodes 
“ Nusselt number 
Static pressure 

- Density weighted probability density f^etion 
~ Pficlet number 

Cell Peclet number 


" Prandtl number 
“ Static pressure 
~ Total pressure 
“ D3mamic pressure 
~ Universal gas constant 
- Streamline radius of curvature 


Reynolds nuinber 

" Generalized source term 
variable ^ ansport equation for 

~ Schmidt number 

- Streamwise. normal, and radial coordinates 
“ Time 

~ Temperature 

- Truncation error in approxinating convection ten. 
Truncation error in approximating diffusion term 
Static temperature 
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Us» Un. 
U 


*i 

Ax, Ay 
Yi 







P 


"‘j 

t 


~ Total temperature 
” Velocity vector 

~ Streamwise, normal, and radial velocity components 

- Radial velocity component 

- Coordinate in i direction (i * 1, 2, 3) 

- Cell dimensions in finite difference forumulation 

- Mass fraction of species i 

~ Chemical production rate of species i per unit mass 

- Exponent of temperature in expression for k.f 

- Kronecker delta 

Turbulent kinetic energy dissipation rate 

- Viscous dissipation correlation in Reynolds stress 
transport equation 

- Eignevalue of 

~ Molecular viscosity 
*• Eigenvalue of Gj 
~ Turbulent eddy viscosity 
~ P + effective eddy viscosity 

- Stoichiometric reactant and product coefficients for 
species i and reaction i 

^ i ■ 

- Mixture fraction (conserved scaler) 

-Mean and variance of the pdf P(5, x) 

- Mixture density 

- Turbulent Pr and tl /Schmidt numbers 

- Drop lifetime 
^ Viscous stress tensor 

- Generalized scalar variable 

- Relaxation factor in iterative 


w 


process 


Ove rbars 


- - Conventionally (time) averaged quantity 

- - Density weighted (time) averaged quantity 

O'* - Conventional fluctuation 

O'*”* “ Density-weighted fluctuation 


HT Nomenc lature 

» Flame emissivity 

* Local average fuel-air ratio 

“ Mean beam length, 1.8 times the dome height, feet 

* Combustor pressure, atmosphere 

■ Gas Temperature, “ R 

■ Luminosity factor, for nonluminous emissivity K * 1.0 
= Film effectiveness 
“ Film Temperature * F 

T coolant * Coolant Temperature * F 

Q/A = Radiation (flame or liner to casing) Btu/hr-ft^ 

^backside “ Cold side heat transfer coefficient Btu/hr-ft^-® p 
hg * Hot side heat transfer coefficient Btu/hr-ft^-“ F 


_f 

a 

Pb 

**8 

Tgas 

K 

n 

**-film 
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